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A.  Problem  Statement 


The  purpose  of  the  proposed  research  was  to  characterize  the  n-th  order 
cyclostationarity  properties  of  general  types  of  LPI  signals  and  to  use  the  characterizations 
to  investigate  and  develop  methods  for  both  secure  communication  and  interception.  For 
secure  communication,  this  includes  means  for  reduction  of  strength  of,  and  elimination  of, 
cyclic  features  that  could  be  exploited  by  an  interceptor.  Signal  design  for  LPI 
communication  was  to  be  considered  from  the  viewpoints  of  both  the  communicator  s 
reception  task  and  the  interceptor's  reception  task.  For  interception,  the  characterizations 
were  to  be  used  to  propose  feature  sets  and  discrimination  rules  for  signal  classification  and 
identification.  The  basic  approach  to  characterizing  n-th  order  cyclostationarity  properties 
of  signals  was  to  be  extended  and  generalized  from  the  recently  developed  theory  of  2nd- 
order  cyclostationarity,  in  which  spectral  characterizations  play  a  crucial  role. 
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B.  Summary  of  Results 


I.  Fundamental  Theory 

In  [1,  2,  3],  the  fundamental  theory  of  higher-order  cyclostationarity  is  presented,  and  some 
applications  of  the  theory  are  discussed.  The  theory  is  an  extension  of  the  theory  of  second- 
order  cyclostationarity,  which  analyzes  the  polyperiodic  components  of  quadratic  transforma¬ 
tions  of  random  time-series,  to  nth-order  cyclostationarity,  which  analyzes  the  polyperiodic 
components  of  nth-order  transformations  of  random  time-series.  The  central  parameters  of 
the  theory  are  ?zth-order  moments  and  cumulants  of  both  the  time-series  itself  (time  domain 
or  temporal)  and  of  finite-time  Fourier  transforms  of  the  time-series  (frequency-domain  or 
spectral).  Unlike  all  other  work  on  cumulants,  in  which  the  cumulant  is  used  because  it  has 
some  well-known  useful  properties,  we  have  shown  that  cumulants  arise  as  the  solution  to 
a  specific  problem  that  is  central  to  the  study  of  higher-order  cyclostationarity.  That  is, 
we  have  shown  that  each  Fourier  component  of  the  ?ith-order  temporal  cumulant  function 
(called  a  cyclic  cumulant)  for  a.  random  polycyclostationary  time-series  can  be  derived  as  the 
solution  to  the  problem  of  calculating  the  part  of  a  sine  wave,  which  is  generated  by  applying 
an  nth-order  homogeneous  polynomial  nonlinearity  to  the  time-series,  that  excludes  contri¬ 
butions  that  are  due  to  the  multiplication  of  sine  waves  that  are  generated  by  the  lower-order 
factors  in  the  homogeneous  polynomial.  In  addition,  Part  I  of  the  papers  [2,  3]  contains  a 
detailed  history  of  the  cumulant  from  its  inception  in  the  late  nineteenth  century  up  to  the 
present.  The  papers  [2,  3]  also  contain  results  on  the  effect  on  the  theory’s  central  parame¬ 
ters  of  certain  signal-processing  operations  that  are  performed  on  the  time-series  in  question, 
including  signal  addition,  signal  multiplication,  convolution,  and  periodic  sampling.  Meth¬ 
ods  for  measurement  of  the  central  parameters  are  presented,  analyzed  mathematically,  and 
evaluated  by  computer  simulation.  The  book  chapter  [4]  presents  a  tutorial  treatment  of 
much  of  the  work  presented  in  [1,  2,  3]. 

II.  Theory  and  Method  for  Measurement  of  Cyclic  Statistics 

In  the  two  conference  papers  [5]  and  [6],  the  measurement  of  the  frequency-domain  cumu- 


2 


lants,  known  as  cyclic  polyspectra,  is  considered.  Much  of  the  material  in  the  first  paper, 
[5],  appears  in  the  later  publications  [1,  2,  3],  but  the  material  in  [6]  has  not  been  pub¬ 
lished  in  any  other  place.  In  [5],  three  methods  of  estimating  the  cyclic  polyspectrum  are 
presented  and  shown  to  converge  to  the  correct  theoretical  functions  when  the  data  record 
length  first  increases  without  bound,  and  then  the  spectral  resolution-width  parameter  is 
allowed  to  approach  zero.  The  correctness  of  the  methods  for  finite  data-record  lengths  and 
spectral  resolution-widths  is  verified  by  computer  simulation.  It  is  also  shown  that  the  three 
measurement  methods  reduce  to  well-known  spectral  estimation  methods  for  order  tw'o  and 
a  cycle  frequency  of  zero.  In  [6],  a  source  of  measurement  error  that  is  particular  to  the  mea¬ 
surement  of  cyclic  polyspectra,  as  opposed  to  polyspectra,  which  are  the  spectral  cumulants 
for  nth-order  stationary  time-series,  is  studied.  This  source  of  error  is  called  leakage  from 
submanifolds,  and  results  from  the  presence  of  lower-order  cyclic  features.  Leakage  from 
submanifolds  is  the  generalization  to  higher  orders  of  the  spectral  leakage  phenomenon  that 
occurs  in  spectrum  estimation  when  the  data  contains  sine-wave  components  (when  the  ideal 
spectrum  contains  impulses).  The  leakage  phenomenon  is  studied  both  analytically  and  by 
computer  simulation  (the  latter  only  for  n  =  4). 

III.  Algorithms  for  TDOA  Estimation 

In  [7,  8],  the  problem  of  estimating  the  time  difference  of  arrival  (TDOA)  between  the 
components  of  a  signal  that  is  received  at  two  spatially  separated  sensors  is  studied  for 
the  case  in  which  the  signal  does  not  exhibit  sufficiently  strong  (exploitable)  second-order 
cyclostationarity  and  is  subject  to  cochannel  interference  and  time- varying  background  noise. 
This  can  happen,  for  example,  when  the  signal  is  a  bandwidth-efficient  communication  signal 
or  when  it  is  a  member  of  a  class  of  low-probability-of-intercept  signals.  In  [7,  8],  work,  a  class 
of  TDOA  estimators  is  derived  to  solve  the  stated  problem.  This  class  of  estimators  is  based 
on  matching  measurements  of  sets  of  nth-order  cross  cyclic  cumulants  in  a  least-squares 
sense.  For  n  =  2,  two  of  the  estimators  in  the  class  reduce  to  known  cyclostationarity- 
exploiting  TDOA  estimators  (SPECCOA  and  CPD).  The  nth-order  version  of  CPD  has 
been  implemented  in  software,  but  simulations  aimed  at  characterizing  its  performance  have 
not  been  completed. 

IV.  A  Method  for  Waveform  Extraction 
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In  [9],  the  theory  of  higher-order  cyclostationarity  is  applied  to  the  problem  of  estimating  a 
signal  waveform  in  the  presence  of  cochannel  interference  and  noise  from  data  obtained  from 
a  single  sensor.  In  particular,  the  signal  of  interest  can  be  completely  spectrally  and  tempo¬ 
rally  overlapped  by  the  interferes.  If  the  signal  exhibits  strong  second-order  cyclostation¬ 
arity,  then  linear  frequency-shift  filtering  can  be  used  to  perform  this  waveform  estimation. 
If  the  interferes  also  exhibit  strong  second-order  cyclostationarity,  then  the  quality  of  the 
waveform  estimate  can  be  improved  by  exploiting  this  cyclostationarity  in  addition  to  that 
of  the  signal  of  interest.  However,  for  the  case  in  which  the  signal  and  interferers  do  not  ex¬ 
hibit  sufficient  second-order  cyclostationarity,  but  do  exhibit  higher-order  cyclostationarity, 
nonlinear  frequency-shift  filtering  can  be  used  to  perform  the  waveform  estimation. 

In  [9],  a  nonlinear  frequency-shift  filter  consisting  of  a  linear  part  and  a  cubic  part  is 
studied  mathematically  and  by  computer  simulation.  The  mathematical  analysis  consists 
of  deriving  the  minimum-mean-squared-error  (MMSE)  design  equations  for  estimating  a 
signal  waveform  contained  in  an  arbitrary  data  record.  The  derived  integral  equations  are 
complicated  and  have  not  yet  been  solved.  The  form  of  the  equations  can  yield  information 
that  helps  to  specify  the  frequency  shifts  to  be  used  in  both  the  linear  and  cubic  parts  of 
the  system.  Such  a  system  was  implemented  in  software  and  simulations  were  performed 
to  assess  the  potential  of  the  method.  In  particular,  it  was  desired  to  determine  if  there 
was  enough  potential  for  solving  problems  of  practical  interest  to  merit  numerical  evaluation 
of  the  MMSE  design  equations.  The  simulations  performed  to  date  show  that,  for  a  small 
number  of  frequency  shifts  and  a  small  number  of  cubic  transformations,  the  reduction  in 
mean-squared  error  (with  respect  to  linear  time-invariant  filtering)  is  about  50%.  However, 
for  the  signal  environments  chosen,  there  is  simply  no  known  viable  alternative  for  waveform 
estimation.  Because  the  method  appears  to  have  some  potential,  it  is  recommended  that 
further  research  be  carried  out  on  nonlinear  frequency-shift  filtering. 

V.  LPI  Signal  Design 

In  [10],  the  challenging  problem  of  designing  a  low-probability-of-intercept  (LPI)  or,  more 
accurately,  a  covert  digital  communication  signal  is  tackled  using  the  theory  of  higher-order 
cyclostationarity.  The  design  philosophy  is  as  follows.  Signals  are  normally  detected  based  on 
their  observable  characteristics,  such  as  power  level,  center  frequency,  bandwidth,  temporal 
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structure,  etc.  When  these  obvious  features  are  hidden  by  a  signal  designer,  as  is  the  case 
with,  for  example,  direct-sequence  spread-spectrum  signals,  other  features  of  the  signal  must 
be  used  to  perform  detection.  A  useful  class  of  features  are  those  periodic  components 
that  can  be  found  in  the  outputs  of  nonlinear  transformations  of  the  signal.  Examples 
of  these  are  symbol  rates,  chip  rates,  and  doubled  carrier  frequencies.  The  interceptor 
nonlinearly  transforms  the  received  signal  and  attempts  to  detect  these  periodic  components 
in  order  to  detect  the  presence  of  the  hidden  signal.  Thus,  the  signal  designer  must  take 
this  interception  strategy  into  account.  He  must  design  his  covert  signal  such  that  no  useful 
periodic  components  can  be  generated  by  a  specified  class  of  nonlinear  transformations  for 
a  given  data-record  length.  Our  LPI  signal  design  starts  at  this  point. 

We  attempt  to  design  a  digital  communication  signal  such  that,  ideally,  no  finite-strength 
additive  sine- wave  components  (with  nonzero  frequency)  can  be  generated  using  nonlineari¬ 
ties  of  order  N  or  less.  (The  ability  to  generate  sine-wave  components  with  zero  frequencies 
cannot  be  avoided,  but  this  does  not  destroy  the  covertness  of  a  signal  either.)  To  minimize 
the  detectability  of  these  sine  waves,  the  digital  communication  signal  is  chosen  to  be  a 
spread-spectrum  signal  and  is,  therefore,  hidden  in  the  much  stronger  noise.  The  theory  of 
second-order  cyclostationary  signals  is  sufficient  to  solve  the  remaining  problem  if  N  is  less 
than  or  equal  to  three,  but  for  larger  Ar,  the  theory  of  higher-order  cyclostationarity  must  be 
used.  This  latter  theory  completely  characterizes  all  the  finite-strength  additive  sine-wave 
components  that  can  be  generated  by  nonlinear  transformation  in  terms  of  cyclic  cumulants. 
We  have  found  that  the  cyclic  cumulants  of  digital  QAM  signals  are  equal  to  the  product 
of  the  cumulant  of  the  symbol  variable  with  a  nonlinear  function  of  the  complex  envelope  of 
the  pulse. 

To  minimize  the  second-order  cyclic  cumulants  of  digital  spread-spectrum  QAM  signals, 
the  chip  envelope  cannot  be  the  standard  rectangular  pulse.  A  zero-percent  excess-bandwidth 
pulse,  such  as  that  used  in  partial-response  signaling,  can  be  used  to  eliminate  all  second- 
order  features  associated  with  the  chip  rate. 

To  minimize  the  second-order  features  associated  with  the  carrier  frequency,  and  all 
higher-order  features,  the  signal  constellation  must  be  appropriately  designed.  This  neces¬ 
sitates  a  departure  from  the  standard  ±1  chip  alphabet.  These  considerations  result  in  a 
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cumulant-minimization  problem  that  involves  a  set  of  nonlinear  equality  constraints  on  the 
chip  alphabet  and  associated  probabilities.  Alternatively,  the  design  problem  can  be  cast  in 
the  form  of  a  large  system  of  nonlinear  equations  resulting  from  setting  all  the  undesirable 
sine-wave  strengths  equal  to  zero,  but  further  research  is  needed  to  determine  if  any  exact 
solutions  to  these  nonlinear  equations  exist. 

In  our  work  to  date,  we  have  developed  numerical  methods  for  solving  the  minimization 
problem  described  above.  These  numerical  methods  have  been  used  to  design  covert/LPI 
signals  with  very  small  second-,  third-,  and  fourth-order  cyclic  cumulants.  Because  of  the 
relationships  between  moments  and  cumulants,  the  corresponding  cyclic  moments  up  to  order 
four  are  also  very  small.  As  an  example,  a  16-QAM  signal  employing  partial-response  pulse¬ 
shaping  has  been  designed  for  which  the  two  largest  sine-wave  strengths  for  nonlinearities 
with  orders  four  and  less  are  17.4  dB  and  28.6  dB  below  the  power  level,  although  the 
strengths  of  most  of  the  regeneratable  sine  waves  are  more  than  30  dB  below  the  power 
level.  By  comparison,  for  BPSK  signals,  the  largest  sine-wave  strength  is  6  dB  above  the 
power  level,  and  for  16-QAM  with  the  conventional  square  constellation  and  rectangular 
pulses,  the  two  largest  are  3.3  dB  and  7.9  dB  below  the  power  level. 

VI.  Cochannel  Signal  Detection 

In  [11],  the  problem  of  determining  the  number  of  cyclostationary  signals  present  (if  any)  in  a 
given  data  record  is  addressed.  The  particular  scenario  of  interest  is  that  for  which  the  signals 
are  temporally  and  spectrally  overlapping.  Lack  of  prior  knowledge,  including  the  number 
of  signals  present,  prohibits  detection  based  on  radiometry  (detection  of  energy  in  certain 
spectral  bands),  and  the  possible  presence  of  signals  that  do  not  exhibit  strong  second-order 
cyclostationarity  rules  out  conventional  cycle  detectors.  However,  an  algorithm  for  blindly 
estimating  cyclic  cumulants  has  been  developed  to  solve  this  problem.  This  algorithm  is 
called  the  general  search  algorithm.  (GSA),  and  its  function  is  to  produce  estimates  of  the 
cyclic  cumulants  and  their  associated  cycle  frequencies  of  the  given  data  for  orders  1  through 
V,  where  N  is  specified  by  the  user  of  the  algorithm. 

The  output  of  the  GSA  consists  of  a  set  of  estimates  of  cyclic  cumulants  and  the  associated 
cycle  frequencies  for  each  order  of  processing.  However,  this  list  of  numbers  is  not  enough  to 
determine  the  number  of  signals  that  are  present  (which  give  rise  to  the  cyclic  cumulants). 
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To  do  this,  the  cyclic  cumulant /cycle  frequency  pair  estimates  need  to  be  grouped  together 
such  that  each  group  corresponds  to  a  single  signal  in  the  data.  The  grouping  algorithm  (GA) 
performs  this  task.  The  GA  has  been  implemented  for  the  case  of  signals  of  the  digital  QAM 
variety,  but  can  be  extended  to  include  classes  of  signals  for  which  formulas  for  the  temporal 
cumulants  are  known.  In  [11],  the  GSA  and  GA  are  described  and  their  performance  is 
illustrated  with  simulations.  This  work,  while  in  an  early  stage  of  research  and  development, 
is  very  promising  and  merits  further  research. 

VII.  Cochannel  Modulation  Identification 

In  [12],  the  GSA  and  GA  algorithms  are  used  to  not  only  detect  any  cyclostationary  signals 
present  in  a  given  data  record,  but  to  identify  their  modulation  types  as  well.  A  catalog  of 
feature  vectors  obtained  by  applying  the  GSA/GA  tandem  to  a  variety  of  noiseless  simulated 
signals  is  obtained.  This  catalog  reveals  that  signals  with  identical  second-order  cyclic  fea¬ 
tures  can  have  radically  different  higher-order  cyclic  features  provided  the  order  N  is  chosen 
large  enough.  For  example,  QPSK  and  8PSK  have  identical  second-order  cyclic  features,  but 
have  distinct  fourth-order  cyclic  features  (QPSK  has  features  related  to  the  carrier  frequency 
for  order  four,  whereas  8PSK  does  not).  To  distinguish  between  8PSK  and  16PSK  requires 
N  =  8.  The  features  are  useful  because  they  can  be  estimated  for  each  signal  in  the  data 
(in  principle)  regardless  of  the  amount  of  spectral  overlap.  An  algorithm  for  identifying  the 
modulation  type  of  a  signal  given  a  measured  feature  has  been  devised  but  has  not  been 
thoroughly  tested.  This  algorithm  measures  the  distance  between  the  measured  feature  and 
each  feature  in  a  set  of  candidate-signal  features  and  chooses  the  modulation  type  that  gives 
rise  to  the  feature  with  minimum  distance. 

VIII.  Nonlinear  System  Identification 

In  [13],  the  problem  of  identifying  the  kernels  of  a  nonlinear  Volterra  system  is  studied.  A 
new  approach  to  obtaining  sets  of  operators  that  are  orthogonal  to  the  Volterra  operators  for 
a  class  of  input  time-series  is  developed.  It  is  shown  how  these  orthogonal  operators  can  be 
used  in  a  crosscorrelation-based  method  to  identify  the  Volterra  kernels  of  nonlinear  time- 
invariant  systems.  This  class  of  inputs  includes  both  Gaussian  and  finite-state  (e.g.,  PSK) 
time-series,  and  simulations  suggest  that  methods  used  with  the  finite-state  inputs  converge 
considerably  more  rapidly  than  the  methods  used  with  Gaussian  inputs,  which  includes  as 
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a  special  case  the  method  of  Wiener,  Lee,  and  Schetzen. 

For  complex-valued  inputs,  the  new  methods  are  believed  to  be  the  only  crosscorrelation- 
based  methods  known  that  identify  arbitrary  order  Volterra  kernels  of  an  infinite-order  sys¬ 
tem  (with  infinite  memory).  Such  methods  can  be  used  to  compute  the  Volterra  kernels  of 
nonlinear  systems  that  can  be  simulated  on  a  computer. 

The  computational  efficiency  of  these  new  crosscorrelation-based  methods  can  be  in¬ 
creased  by  using  an  FFT  algorithm  to  implement  the  frequency-domain  counterparts  devel¬ 
oped  in  this  project  that  are  based  on  frequency-smoothed  cross-spectra. 

IX.  Polyperiodic  Nonlinear  System  Identification 

In  [14],  the  class  of  polyperiodic  nonlinear  (PPN)  systems  that  admit  a  new  double  series 
representation,  called  the  Fourier- Volterra  series,  is  considered,  and  a  class  of  cyclostationary 
random  inputs  that  enables  the  analytical  specification  of  sets  of  operators  on  the  input  that 
are  orthonormal  over  all  time  to  the  Fourier- Volterra  operators  are  defined.  The  reciprocal 
operators  are  used  to  obtain  a  crosscorrelation  formula  for  identifying  the  Fourier- Volterra 
kernels  of  PPN  systems.  It  is  shown  how  the  orthonormal  operators  for  the  finite-state 
phase-modulated  and  Gaussian  amplitude-modulated  inputs  can  be  determined.  The  finite- 
state  inputs  are  expected  to  result  in  computationally  attractive  identification  relative  to 
the  computational  load  required  for  Gaussian  inputs  because  of  the  substantially  longer 
averaging  times  that  are  apparently  required  by  the  latter.  Nevertheless,  the  computational 
load  for  the  former  can  be  very  substantial  not  only  because  of  the  averaging  time  required, 
but  also  because  of  the  number  of  terms  required  in  the  double  series  representation.  Some 
reduction  in  the  computational  load  of  the  crosscorrelation  method  can  be  achieved  by 
using  an  FFT  algorithm  to  implement  its  frequency-domain  counterpart,  which  is  based  on 
frequency-smoothed  cyclic  cross-periodograms,  and  which  was  developed  in  this  project. 
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Design  of  Low-Probability-of-Intercept  Signals: 
Minimum-Cumulant  Digital  QAM  Signals 


Chad  M.  Spooner  and  William  A.  Gardner 
July  8,  1994 

Abstract 

Effective  methods  of  signal  interception  for  signals  that  are  hidden  in  noise  include 
those  that  are  based  on  exploiting  second-order  cyclostationarity.  The  exploitability 
of  second-order  cyclostationarity  for  interception  purposes  can  be  measured  by  the 
peak  strength  of  the  sine-wave  components  that  appear  when  the  signal  is  subjected  to 
certain  quadratic  nonlinearities.  To  counter  such  cyclostationarity-based  interception, 
signals  can  be  designed  to  yield  quadratically  regeneratable  sine  waves  with  small 
strengths.  Thus,  interceptors  are  forced  to  use  other  means  to  intercept  such  signals. 
One  class  of  interception  methods  for  these  LPI  signals  is  based  on  the  exploitability 
of  higher-order  (order  larger  than  two)  cyclostationarity.  For  this  interception  to  be 
effective,  the  peak  strengths  of  the  sine-wave  components  in  the  output  of  higher-order 
nonlinearities  must  be  large  enough  to  detect.  The  signal  designers  are  then  forced  to 
minimize  these  strengths.  In  this  work,  digital  QAM  signals  with  minimum  second-, 
third-,  and  fourth-order  cyclostationarity  are  designed  using  a  combination  of  pulse¬ 
shaping  and  amplitude/phase  constellation  design  based  on  the  newly  developed  theory 
of  higher-order  cyclostationarity. 


1  Introduction 

A  substantial  amount  of  recent  research  work  on  signal  interception  has  focussed  on  exploit¬ 
ing  the  cyclostationary  nature  of  digital  communication  signals,  such  as  pulse-amplitude 
modulation  (PAM),  phase-shift  keying  (PSK),  quadrature-amplitude  modulation  (QAM), 
direct-sequence  spread-spectrum  signals,  and,  to  a  lesser  degree,  frequency  modulation, 
frequency-shift  keying,  and  frequency-hopped  spread-spectrum  signals  [5]— [9] .  The  key  to 
such  interception  is  the  ability  to  generate  a  sine  wave  (or  set  of  sine  waves)  by  nonlinearly 
transforming  the  signal,  and  then  detecting  the  presence  of  this  regenerated  sine  wave.  This 
interception  can  be  accomplished  even  when  the  background  noise  and  interference  is  very 
strong  provided  that  enough  data  is  available  [1],  Thus,  the  degree  of  detectability  of  a  cyclo¬ 
stationary  signal  can  be  measured  in  terms  of  the  strength  of  its  cyclostationarity,  which  can 
be  defined  in  terms  of  the  strengths  of  its  regeneratable  sine  waves  for  a  fixed  signal-power 
level  [6]. 

These  interception  considerations  have  necessarily  had  an  impact  on  the  design  of  low- 
probability-of-intercept  (LPI)  signals.  In  particular,  it  is  desirable  to  limit  the  amount  of 
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cyclostationarity  that  can  be  exploited  by  the  interceptor,  while  simultaneously  allowing 
the  communicator  to  detect  and  demodulate  the  signal  by  using  a  receiver  with  reasonable 
complexity.  Thus,  it  is  desired  to  minimize  the  strengths  of  the  regeneratable  sine  waves, 
but  because  regenerated  sine  waves  are  often  used  (explicitly  or  implicitly)  by  the  intended 
receiver  to  synchronize  to  the  signal,  this  is  not  just  a  simple  minimization  problem;  there  is 
a  tradeoff  between  probability  of  intercept  and  difficulty  of  reception  for  the  intended  com¬ 
municator.  Nevertheless,  in  this  report  the  focus  is  on  the  minimization  problem  without 
regard  to  the  resulting  complications  for  the  intended  communicator.  The  idea  is  to  deter¬ 
mine  the  limits  on  wffiat  can  be  done  by  using  the  proposed  approach,  and  then  to  evaluate 
the  difficulty  of  reception  of  the  designed  signals. 

One  way  to  minimize  the  strength  of  regeneratable  sine  waves  that  applies  to  the  class 
of  digital  QAM  signals  is  through  pulse-shaping.  That  is,  since  digital  QAM  signals  (at 
radio  frequency)  can  be  represented  by  complex-valued  PAM  signals  (at  baseband),  then 
their  spectral  correlation  functions  are  easily  characterized  in  terms  of  the  baseband  pulse 
transform  [1].  The  support  of  this  transform  determines  the  number  of  regeneratable  sine 
waves  as  well  as  their  strengths.  By  bandlimiting  the  pulse  transform  to  the  Nyquist  band 
(where  the  positive-frequency  bandwidth  of  the  pulse  is  equal  to  half  the  pulse  rate),  the 
second-order  cyclostationarity  of  the  signal  that  is  associated  with  the  pulse  rate  is  annihi¬ 
lated  [1].  The  radio  frequency  signal  can  still  exhibit  cyclostationarity  associated  with  the 
carrier  if  the  symbol  constellation  (discrete  support  of  the  probability  mass  function  of  the 
complex-valued  random  symbols)  does  not  meet  certain  symmetry  requirements. 

Pulse-shaping  can  be  used  to  reduce  the  exploitable  second-order  cyclostationarity  of  a 
signal,  but  it  is  less  effective  at  reducing  the  exploitable  higher-order  cyclostationarity  of  the 
signal.  The  positive-frequency  bandwidth  of  the  pulse  transform  would  need  to  be  reduced  to 
1/n  times  the  pulse  rate  (or,  more  generally,  the  support  of  the  pulse  transform  would  have  to 
be  limited  in  some  other  way)  in  order  to  annihilate  the  ?rth-order  regeneratable  sine  waves, 
which — for  n  =  2 — would  introduce  an  unacceptable  amount  of  intersymbol  interference  to 
the  transmitted  signal,  thus  foiling  the  intended  receiver.  However,  the  symbol  constellation 
itself  can  be  designed  to  minimize  the  strength  of  the  nth-order  regeneratable  sine  waves. 
This  is  a  result  of  the  fact  that  the  strengths  of  the  pure  nth-order  sine  waves  [4]  for  a 
given  digital  QAM  signal  are  directly  proportional  to  the  nth-order  cumulant  of  the  symbol 
constellation. 

This  report  discusses  the  problem  of  minimizing  the  nth-order  cumulants  of  discrete 
random  variables  by  designing  their  probability  mass  functions  for  the  purpose  of  design¬ 
ing  signal  constellations  and  associated  probabilities  for  a  digital  QAM  (DQAM)  signaling 
scheme.  The  cumulant-minimization  problem  is  also  interesting  from  a  mathematical  view¬ 
point  because  it  can  be  thought  of  as  an  attempt  to  find  a  small-alphabet  discrete  random 
variable  that  has  a  PDF  that  is,  in  some  suitable  sense,  a  good  approximation  to  the  con¬ 
tinuous  Gaussian  probability  density  function  (PDF).  This  follows  because  the  higher-order 
cumulants  for  Gaussian  random  variables  are  all  zero.  In  fact,  this  observation  is  the  basis 
for  one  of  the  numerical  methods  described  in  this  report. 

In  the  following,  the  cumulants  (and  moments)  of  digital  QAM  signals  are  presented,  and 
then  the  minimization  problem  is  stated  and  explored  analytically.  Numerical  approaches  to 
this  minimization  problem  are  described,  and  some  results  are  provided  and  discussed.  The 
reader  of  this  report  is  assumed  to  be  familiar  with  higher-order  statistics  in  general,  and 
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with  higher-order  cyclostationarity  in  particular  [4]. 


2  Moments  and  Cumulants  of  PAM  Signals 


The  class  of  signals  considered  herein  is  specified  by 
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where  3?{-}  is  the  real-part  operator,  fc  is  the  carrier  frequency,  9C  is  the  carrier  phase,  T0 
is  the  pulse  interval  length,  t0  is  the  pulse  timing  constant,  and  cm  and  sm  are  real-valued 
independent  symbol  variables  that  take  on  values  in  some  finite  set.  The  analytic-signal 
representation  of  s(t)  is  denoted  by  sa(t),  and  is  given  by 


sa(t)  =  E 


771— —  CX>  L 


e  c 


2  (^ti  ISm ) 


p(t  +  mT0  +  t0)e 


i2irjct 
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and  the  complex  envelope  of  the  DQAM  signal  s(t)  is  denoted  by  sc(t): 

J6c 
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oo 


p(t  +  mT0  +  t0), 


=  E  a,mp(t +  mT0  +  t0), 


(3) 


which  is  a  complex-valued  PAM  signal.  This  signal  model  includes  PSK,  M ary-QAM,  and 
(more  generally)  amplitude-and-phase-shift-keying  (APK)  modulations. 

The  cumulants  of  a  complex  PAM  signal  sc(t)  with  independent  and  identically  dis¬ 
tributed  symbols  are  given  by  [4] 


CSc(t,  7_)n 


Cumulant 
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E 


m——oo 


fl  +  ™T0  +  Tj  + 10)  , 
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where  (*)j  denotes  optional  conjugation  of  the  j th  element  in  the  set  {si*^(t  +  i,  and 

Ca,n  is  the  nth-order  cumulant  of  the  symbol  sequence  {am},  and  is  given  by 


in  which 
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and  Pn  is  the  set  of  all  distinct  partitions  \v i,  •  *  - ,  np}  of  the  set  of  indices  {1, 2,  •  -  • ,  n},  and  p 
is  the  size  of  the  partition  (1  <  p  <  n)  [2].  The  reduced-dimension  cyclic  temporal  cumulant 
function  (RD-CTCF)  and  cyclic  polyspectrum  (CP)  follow  directly  from  (4): 


?i(fX  =  0  =  k/TOl  (8) 

4o  j—  i 


where  P(f)  is  the  Fourier  transform  of  p(t).  Thus,  the  nth-order  temporal  and  spectral 
cumulants  of  the  PAM  signal  are  directly  proportional  to  the  nth-order  cumulant  of  the 
symbol  variable. 

The  cyclic  temporal  moment  functions  (CTMFs)  of  the  PAM  signal  sc(t)  can  be  found 
by  combining  the  cumulants  of  sc(t)  in  the  following  alternative  ways  [4]: 


RZ(r)n  =  C“(r)„  +  Y. 


Pn 

p^l 


Pn 
p¥=  1 


L  /3*W=1 


(-i)p_i(p  - 1)!  £  n«2K)»i 

athoh1 
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where  a  and  (3  are  vectors  of  impure  and  pure  cycle  frequencies,  respectively  [4].  These 
expressions  are  complicated  for  n  >  3,  and  no  simplifications  have  been  found.  Clearly,  the 
cyclic  moment  is  not  generally  easy  to  compute,  and  it  is  obviously  not  proportional  to  the 
moment  of  the  symbol  variable  Ra:Tl  in  general.  Nevertheless,  for  the  case  in  which  p(-)  has 
width  less  than  or  equal  to  To,  and  Tj  =  u  for  all  j,  the  CTMF  is  given  by 


/?  r°°  n  /x 

<(«)»  =  ^  /  II  PW'(<  +  u)e~'2M dt  e'2nato, 

1  0  J-oo  j_  J 

for  a  =  k/To-  Unfortunately,  for  rectangular  pulses  this  particular  CTMF  is  zero  for  k  ^  0 
(but  it  is  nonzero  for  other  r).  In  general,  then,  it  is  difficult  to  relate  the  CTMF,  which 
is  the  strength  of  a  sine-wave  component  of  an  nth-order  lag  product,  to  the  probabilistic 
parameters  of  the  signal  constellation  for  arbitrary  n  and  r,  whereas  it  is  easy  to  relate 
the  CTCF,  which  is  the  strength  of  a  pure  nth-order  sine  wave  component  in  an  nth-order 
lag  product,  to  these  parameters.  However,  if  there  is  no  lower-order  cyclostationarity, 
then  R°c (r)n  =  C"(r)„,  and  the  strength  of  the  pure  nth-order  sine  wave  is  equal  to  the 
strength  of  the  impure  nth-order  sine  wave.  Thus,  in  this  latter  case,  both  strengths  are 
directly  proportional  to  Ca,n ■  The  most  common  example  of  this  is  the  case  of  n  =  4  and 
0%  excess-bandwidth  pulses  (for  example,  partial  response  signals  [3],  of  which  duobinary 
signals  are  an  example).  For  higher  orders,  it  is  much  less  likely  that  there  will  be  no  lower- 
order  cyclostationarity  and,  therefore,  the  impure  nth-order  sine  waves  are  more  difficult  to 
analyze.  Nevertheless,  we  shall  discover  that  designing  a  signal  to  exhibit  minimum  second- 
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and  fourth-order  cyclic  cumulants  results  in  a.  signal  that  has  small  second-  and  fourth-order 
cyclic  moments  as  well. 

In  general,  cyclic  cumulants  are  signal  selective  and  tolerant  to  Gaussian  corruption, 
whereas  cyclic  moments  are  not  [2,  4].  In  addition,  by  inspection  of  the  tables  to  be  pre¬ 
sented,  the  cumulants  of  the  signal  constellation  get  very  large  with  increasing  n  whereas  the 
moments  do  not  (for  a  fixed  constellation  variance).  Thus,  an  interceptor  might  be  wise  to 
compute  cumulants  instead  of  moments,  and  an  LPI  signal  designer  might  be  wise  to  focus 
on  minimizing  the  exploitability  of  cumulants  rather  than  moments.  All  these  considerations 
lead  us  to  focus  on  the  cumulants  of  the  signal  constellation,  rather  than  the  moments. 


3  Moments  and  Cumulants  of  Symbol  Constellations 

To  better  understand  the  nature  of  the  cumulants  and  moments  of  the  complex  PAM  signal 
sc(t),  the  values  of  Ca<n  and  Ra^n  corresponding  to  some  simple  DQAM  signals  are  presented 
in  this  section.  The  values  of  am  (the  signal  constellation  and  associated  probabilities)  for 
0C  =  0  for  BPSK,  4-level  PAM,  QPSK,  8PSK,  16PSK,  32PSK,  64PSK,  4QAM,  8QAM, 
16QAM,  and  64QAM  signal  types  are  tabulated  in  Table  1.  The  constellations  are  scaled  so 
that  their  variance  is  equal  to  one. 

For  a  few  distributions  of  am  (constellations),  the  nth-order  cumulants  can  be  computed 
analytically,  but  the  resulting  expressions  always  contain  a  number  (a  function  of  n)  that  can 
only  be  computed  by  means  of  a  recursion  or  an  infinite  series.  The  values  of  moments  and 
cumulants  given  in  the  tables  herein  were  computed  using  a  computer  program  written  in  C, 
and  are  shown  in  Tables  2,  3,  and  4.  This  program  can  compute  the  nth-order  moments  and 
cumulants  of  any  finite-alphabet  random  variable.  The  results  for  the  analytically  tractable 
cases  match  the  corresponding  results  in  the  tables.  Several  analytical  results  concerning  the 
cumulants  of  discrete  random  variables  are  given  in  Section  5.  These  moments  and  cumulants 
also  indicate  the  relative  strengths  of  the  regeneratable  sine  waves  for  the  various  signal  types 
that  are  listed.  For  example,  the  cumulants  of  BPSK  grow  most  rapidly  with  the  order  n, 
indicating  that  this  signal  is  the  most  detectable  of  those  listed  in  the  tables:  in  some  sense 
the  ±1  binary  random  variable  is  the  opposite  of  the  Gaussian  random  variable.  We  have 
found  no  other  random  variable  with  larger  cumulants.  This  observation  is  consistent  with 
the  fact  that  BPSK  signal  also  exhibits  the  strongest  second-order  cyclostationarity  of  any 
known  modulated  signal. 


4  Cycle  Frequencies  of  Digital  QAM  Signals 

From  (7)  and  (8),  the  cycle  frequencies  of  the  complex-valued  PAM  signal  (the  complex 
envelope  of  the  DQAM  signal)  are  restricted  to  harmonics  of  the  pulse  rate  1/T0,  and  the 
width  of  the  Fourier  transform  of  the  pulse  function  p(t)  determines  the  number  of  such 
harmonics  that  result  in  nonzero  values  of  Cfc(-)„.  We’ll  not  concern  ourselves  with  the 
latter  detail  for  now.  We  simply  find  the  cycle  frequencies  for  the  analytic  signal  assuming 
that  the  complex  envelope  exhibits  cyclostationarity  at  all  harmonics  of  the  pulse  rate.  This 
is  the  case,  for  example,  when  the  function  p(t)  is  a  rectangle  with  width  T0  because  in  this 
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case  the  Fourier  transform  P(f)  has  infinite  width.  These  results  can  be  easily  modified  if 
the  width  (support)  of  P(-)  is  finite:  there  will  be  an  upper  limit  to  k  in  |fc|/To  beyond  which 
all  cyclic  cumulants  are  zero. 

The  cycle  frequencies  for  the  complex  envelope  for  a  particular  modulation  type  can  be 
determined  by  using  the  tables  of  values  for  the  cumulants  Ca,n-  When  Ca,n  =  0,  there  are  no 
cycle  frequencies  (because  there  are  no  features),  and  when  Ca,n  /  0,  the  cycle  frequencies 
are  simply  the  harmonics  of  the  pulse  rate. 

Let  A™  denote  the  cycle  frequencies  of  the  complex  envelope  sc(t)  for  order  n  with  m 
factors  conjugated.  It  can  be  shown  by  using  results  in  [2]  that  the  cycle  frequency  set  B '™ 
for  the  analytic  signal  is  given  by 


B 


m 

n 


{a  :  a 


7  +  (n  -  2m) fe,  7  €  A™}  , 

0, 


A™  /  0, 

otherwise. 


The  cycle  frequencies  for  the  analytic-signal  representations  of  some  common  digital  com¬ 
munication  signals  are  shown  in  Table  5. 


5  Cumulants  of  Certain  Discrete  Random  Variables 

The  cumulants  of  certain  discrete  random  variables  can  be  obtained  by  using  the  definition  of 
the  cumulant  as  the  derivative  of  the  logarithm  of  the  characteristic  function  of  the  variable, 
say,  A: 

Ca,«  =(*')■”  In  • 

Consider  a  binary  symmetric  real-valued  random  variable  A  with  PDF 

Ja{u )  =  +  a)  +  5(u  -  a)]. 

The  characteristic  function  for  A  is 


$4(0;)  =  cos(urn), 

and  the  nth-order  cumulant  is  given  by 
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where  Bm  is  the  mth  Bernoulli  number, 
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the  nth-order  cumulant  for  A  is  zero  for  odd  n,  and  for  even  n  it  simplifies  to 

an2n( 2”  -  1  )Bn 


CA,n  = 


n 


(12) 


Next,  consider  a  quaternary  symmetric  variable  with  density  function 


fA{u)  =  ^-[<S(u  +  a)  +  6(u  -  a)  +  5(u  +  b)  +  £(u  -  £>)]. 


The  characteristic  function  for  A  is 

$4(0;)  =  ^[cos(tua)  +  cos(u;6)]. 

By  using  a  trigonometric  identity,  the  log-characteristic  function  can  be  expressed  as 
In  $4(0;)  =  In  cos(cu[a  +  b\/2)  +  lncos(u;[a  -  b]/ 2). 

It  follows  that  the  cumulant  for  A  is  given  by 

CU,„  =  [(«  +  6)”  +  (a  -  (>)”)  (13) 

for  even  values  of  n,  and  is  equal  to  zero  for  odd  values  of  n.  By  looking  up  the  values  of 
the  Bernoulli  numbers,  the  values  listed  in  Table  2  (for  BPSK  and  4-level  PAM),  which  are 
obtained  by  numerical  evaluation,  can  be  confirmed1. 


6  Minimizing  Cumulants 

The  cumulants  of  a  binary  symmetric  random  variable  (12)  are  nonnegative  and  therefore 
their  minima  must  correspond  to  the  symbol  values  of  a  =  —  a  =  0  (cf.  (12)).  Similarly, 
for  the  quaternary  symmetric  variable,  the  minimum  cumulant  (13)  for  each  n  corresponds 
to  a  =  b  =  0.  However,  if  we  incorporate  the  constraint  that  the  variance  of  the  variable 
must  equal  some  constant,  say  1,  then  it  is  not  clear  that  there  is  no  useful  minimum.  For 
example,  for  the  symmetric  quaternary  random  variable,  if  the  variance  is  equal  to  1,  then 
the  extrema  of  the  fourth-order  cumulant  are  given  by  the  solutions  to  the  following  two 
equations 

(a  ±  \/2-a2)3(l  T  )  +  (a  =F  ^2  -  a2)3(l  ±  ^==)  =  0. 

These  equations  are  satisfied  when  a  =  ±1,  and  a  =  0.  The  former  two  solutions  result 
in  the  random  variables  a\  =  1,  b\  =  —  1  with  probabilities  of  0.5,  and  the  latter  solution 

xThe  first  several  Bernoulli  numbers  are  given  by  B2  =  1/6,  B4  =  —1/30,  Be  =  1/42,  B$  =  —1/30. 
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results  in  the  random  variable  eq  =0,6:  —  \/2 ,b2  =  —  \/2  with  probabilities  0.5,  0.25,  and 
0.25,  respectively.  The  fourth-order  cumulant  of  the  former  random  variable  is  given  by  (12) 
with  n  =  4  and  a  =  1  (which  equals  —2),  and  that  for  the  latter  random  variable  is  given  by 

2\/24 B4(24  -  1)  _  x 

4 

These  two  random  variables  are  not  particularly  useful  because  the  associated  DQAM  signals 
still  have  relatively  strong  fourth-order  cyclostationarity.  It  follows,  then,  that  to  obtain 
meaningful  results  from  a  minimization  procedure,  we  must  remove  some  of  our  constraints 
on  the  random  variable  distribution,  that  is,  we  must  consider  asymmetric  random  variables 
and/or  those  with  unequal  probabilities.  In  the  case  of  the  quaternary  symmetric  random 
variable,  we  might  consider  instead  a  quaternary  equiprobable  random  variable  with  four 
arbitrary  values,  or  we  can  allow  the  probabilities  to  deviate  from  0.25,  or  both.  In  these 
cases,  the  resulting  formulas  for  the  cumulants  are  difficult  to  analyze,  and  we  must  resort 
to  numerical  methods. 

In  general,  for  a  constellation  with  M  real-valued  symbols,  the  formula  for  Ca,n  is  an 
nth-order  polynomial  in  2 M  variables  (although  one  of  the  probability  variables  can  be 
eliminated).  Therefore,  to  minimize  \Ca,n\  requires  the  simultaneous  solution  of  2 M  poly¬ 
nomial  equations,  with  maximum  order  n.  As  seen  from  the  example  above,  additional 
constraints  (such  as  unit  variance)  must  be  incorporated  for  the  results  to  be  meaningful. 
Because  of  the  large  number  of  variables  and  the  nonlinear  nature  of  the  equations,  this 
problem  appears  to  be  generally  too  difficult  to  solve  analytically,  and  so  numerical  methods 
must  be  used. 


7  Constraints  for  Complex  Symbols 

In  general,  we  consider  the  problem  of  finding  a  four-symbol  signal  constellation  with  mini¬ 
mum  the  fourth-order  cumulant  subject  to  the  following  constraints  on  the  symbol  variable: 

1.  Zero  mean 

2.  Unit  variance 

3.  No  two  symbols  are  the  same  (unique  symbols) 

4.  No  symbol  has  zero  probability 

5.  Certain  second-,  third-,  and  fourth-order  moments  are  zero. 

The  signal-design  considerations  that  lead  to  these  constraints  are  explained  in  detail  next. 

In  order  for  the  designed  DQAM  signal  to  be  useful,  the  mean  of  the  constellation  must 
be  zero  for  power  efficiency.  To  compare  the  cumulants  of  the  candidate  constellations,  it 
is  necessary  to  normalize  them  somehow.  The  normalization  chosen  is  unit  variance.  This 
constraint  is  important  because  the  cumulants  can  be  made  very  small  (i.e.,  zero)  by  setting 
all  the  symbols  to  zero,  and  this  trivial  solution  must  be  avoided.  In  order  for  the  designed 
DQAM  signal  to  be  useful,  the  symbols  must  be  unique,  no  symbol  can  have  zero  probability, 
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and  the  symbol  magnitude  cannot  exceed  some  constant.  If  a  symbol  has  zero  probability, 
the  constellation  has  fewer  than  the  desired  number  of  elements,  that  is,  the  bits  can  be 
coded  with  fewer  than  the  specified  number  of  symbols.  The  problem  of  finding  a  minimum 
cumulant  for  this  smaller  constellation  size  is  deemed  a  distinct  problem.  If  the  symbols 
are  not  unique,  then  again  we  have  a  smaller-size  constellation.  Finally,  the  bound  on  the 
symbol  magnitude  is  necessary  to  bound  the  dynamic  range  of  the  resultant  DQAM  signal. 

Let  denote  the  complex- valued  symbols  of  the  complex-envelope  signal  repre¬ 

sentation  of  the  passband  signal.  Symbol  aj  occurs  with  probability  pj.  The  foregoing 
constraints  can  be  expressed  mathematically  as  follows: 


M 

ajPj  —  0) 

(14) 

j=i 

M 

E  M  2Pj  =  h 

(15) 

j= 1 

0  <  \a3\  <  K,  1  <3  <  M, 

(16) 

1  <  Pj  <  £2,  1  <j  <  M. 

(17) 

The  remaining  constraints  follow  from  forcing  the  carrier-related  cyclic  cumulant  features 
for  orders  two,  three,  and  four  to  zero,  and  then  finding  the  constellation  that  minimizes  the 
fourth-order  symbol-rate  features  by  minimizing  the  appropriate  fourth-order  cumulant. 

The  second-order  cyclostationarity  (SOCS)  that  is  associated  with  the  symbol-rate  cycle 
frequencies  vanishes  because  the  pulse  is  assumed  to  have  positive-frequency  bandwidth  less 
than  or  equal  to  half  the  pulse  rate.  The  SOCS  that  is  associated  with  the  carrier  frequency 
vanishes  if  the  second-order  cumulant  of  the  symbol  variable  with  no  conjugations  is  zero, 
which  leads  to  the  fifth  constraint: 

M 

Ca,2  =  £«?W  =  °.  (18) 

j=  1 

Equation  (18)  follows  from  (14);  that  is,  the  second-order  cumulant  is  equal  to  the  second- 
order  moment  if  the  mean  is  zero.  Constraint  (18)  implies  that  the  second-order  cumulant 
with  both  factors  conjugated  (C„  2)  is  also  zero.  Notice  that  the  notation  C%p  is  used  to 
denote  the  pth-order  cumulant  of  a  with  q  conjugated  factors. 

The  third-order  features  are  forced  to  zero  in  a  similar  manner.  Because  the  third-order 
cumulants  of  the  symbol  variable  are  equal  to  the  third-order  moments  with  certain  products 
of  lower-order  moments  subtracted  off  (each  of  which  contains  at  least  one  occurrence  of  the 
mean),  the  two  constraints  on  the  third-order  cumulants  reduce  to  the  following: 

M 

C”  3  =  (cy  =  £  “h  =  0,  (19) 

j= i 

M 

Cl3=(Cf,3r  =  £«Kft  =  0.  (20) 

j= 1 
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Similarly,  the  carrier-related  fourth-order  cyclic  cumulant  features,  which  occur  for  0,  1, 
3,  and  4  conjugated  factors  are  forced  to  zero  by  the  following  two  constraints: 

M 

C  =  (O' =  Edft  =  °.  (21) 

j= i 

M 

O  =  (0‘  =  £«?“&  =  «•  (22) 

.)  =  ! 

If  all  these  constraints  are  met,  then  the  remaining  cyclostationarity  of  order  four  or  less  can 
only  occur  for  order  four  with  two  conjugated  factors,  which  is  cyclostationarity  associated 
with  the  pulse  rate.  The  goal  of  the  numerical  processing  is  to  find  the  symbols  and  associated 
probabilities  that  meet  the  constraints  ( 14)— (22)  and  minimize  the  fourth-order  cumulant 
that  is  associated  with  the  symbol-rate  cyclic  features: 

M 

cIa  =  E  Kl  ~  2 

3= 1 

8  Constraints  for  Real  Symbols 

As  with  complex  symbols,  the  constraints  on  real  symbols  can  be  determined  by  remembering 
that  the  symbols  correspond  to  a  PAM  signal  that  multiplies  a  complex-valued  carrier  wave. 
It  is  desired  to  eliminate  the  cyclostationarity  of  the  modulated  signal.  This  can  be  done  by 
minimizing  the  first-,  third-,  and  fourth-order  cumulants  of  the  PAM  signal  by  minimizing 
the  cumulants  of  the  symbol  variable.  Constraints  (14)— (17)  apply  unchanged.  The  SOCS 
associated  with  the  doubled-carrier  frequency  cannot  be  eliminated  from  the  signal  without 
forcing  the  power  or  bandwidth  to  zero.  The  third-order  cyclostationarity  is  destroyed  by 
constraint  (19).  The  objective  is  to  minimize  the  fourth-order  cumulant  of  the  symbol 
variable: 

M 

CaA  =  E  ~  3 
j=l 

subject  to  ( 14)— (17)  and  (19). 

9  Digital  QAM  Signals  From  Real  Random  Variables 

Although  the  real-valued  distributions  obtained  by  meeting  the  constraints  described  in  the 
previous  section  cannot  directly  be  used  to  construct  a  complex-valued  digital  QAM  signal 
with  small  second-,  third-,  and  fourth-order  cyclic  features,  they  can  be  used  indirectly  to 
do  so,  that  is,  the  resulting  signal  has  a  complex-valued  PAM  representation  that  meets  the 
constraints  of  Section  7,  with  the  single  exception  of  (21).  Nevertheless,  the  cumulant  C°4 
is  made  as  small  as  the  cumulant  that  is  minimized  (C%  2),  as  explained  next. 

Let  b  and  c  be  independent  and  identically  distributed  real-valued  random  variables  that 
meet  the  constraints  of  the  previous  section,  and  let  these  random  variables  have  fourth-order 


M 

E  ajPJ 

3  =  1 


M 

£ 

3= 1 


«jl  Pj 


(23) 
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cumulant  C .  The  results  of  Appendix  A  imply  that  the  complex-valued  random  variable 
a  —  (b  T  ic)/y/2  has  the  following  cumulants: 


/>0 

°a,l 

=  0 

°a,2 

=  0 

C1 

'-'a, 2 

=  1 

/i0 

Ua,3 

=  0 

Cl* 

=  0 

Cl  4 

=  C/2 

G1 

Ua,4 

=  0 

r 2 

L/a,4 

=  <7/2. 

It  follows  that  minimizing  the  fourth-order  cumulant  of  b  minimizes  the  fourth-order  cumu¬ 
lants  of  a. 


10  Numerical  Methods  of  Minimization 

The  two  numerical  methods  used  to  find  random  variables  with  minimum  fourth-order  cum¬ 
ulants  are  described  in  this  section. 

10.1  Exhaustive  Search 

This  program  searches  for  the  real-  or  complex-valued  symbol  constellation  that  gives  the 
smallest  nth-order  cumulant  by  exhaustive  search.  The  range  of  symbol  values  and  the 
fineness  of  the  search  grid  are  among  the  input  parameters.  An  optional  input  parameter 
is  an  initial  constellation,  so  that  results  from  a  coarse  search  can  be  refined.  The  symbol 
variable  range  and  the  probabilities  are  discretized  and  the  constraints  are  checked  for  each 
possible  combination  of  symbols  and  probabilities  (with  the  additional  constraint  that  the 
probabilities  must  sum  to  one).  If  the  constraints  are  met,  then  the  fourth-order  cumulant  is 
computed  and  compared  to  the  minimum  obtained  thus  far.  This  program  is  computationally 
expensive  for  random  variables  with  large  alphabet  size,  and  can  only  be  used  to  find  small- 
alphabet  distributions  by  using  the  currently  available  computing  resources. 

10.2  Gaussian  Approximation 

The  real-valued  discrete  random  variable  with  uniform  spacing  and  with  minimum  nth-order 
cumulant  obtained  by  approximating  the  Gaussian  PDF  is  sought.  By  uniform  spacing  we 
mean  that  the  difference  between  every  two  adjacent  symbols  (on  the  real  line)  is  a  constant. 
This  approximation  is  specified  by  two  parameters:  width  and  size.  The  width  parameter  is 
the  width  of  the  Gaussian  distribution  that  is  split  into  size  pieces,  the  centers  of  which  are 
uniformly  spaced.  Thus,  the  first  element  in  the  constellation  corresponds  to  the  amount 
of  area  in  the  tail  of  the  Gaussian,  and  subsequent  elements  have  probability  equal  to  the 
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area  of  the  Gaussian  centered  at  the  element  value,  with  width  equal  to  width/ size.  The 
parameter  size  is  equal  to  the  size  of  the  desired  signal  constellation. 

This  approximation  method  is  an  attempt  to  answer  the  question:  Where  does  one  place 
size  impulses  on  the  real  line  such  that  the  higher-order  cumulants  for  the  resultant  PDF 
are  minimum? 

11  Numerical  Results 

11.1  Real- Valued  Constellations 

In  the  case  of  real-valued  variables,  the  exhaustive  search  method  is  used  for  the  case  of 
M  —  4,  and  the  Gaussian  approximation  method  is  used  for  M  =  4,8,  and  16. 

In  the  case  of  exhaustive  search,  the  four  symbols  are  constrained  to  have  magnitude 
greater  than  or  equal  to  0.0,  but  less  than  or  equal  to  I\  —  3.0.  For  a  symbol  grid  with  fineness 
6a  and  probability  grid  with  fineness  Sp  =  e 2,  the  results  of  the  computer  searches  using  the 
exhaustive  search  method  are  shown  in  Table  6.  For  the  case  of  equiprobable  symbols, 
8a  =  0.05,  and  the  same  bounds  on  the  symbol  magnitude,  the  resulting  distribution  is 
±1.4,  ±0.15,  and  the  cumulant  corresponding  to  this  distribution  has  magnitude  1.03.  Note 
that  this  distribution  is  close  to  that  obtained  by  analysis  in  Section  6. 

The  Gaussian  approximation  technique  yields  the  minimum  cumulants  and  corresponding 
constellations  shown  in  Table  7.  For  the  case  of  M  —  4,  preliminary  searches  showed  that 
the  optimal  width  parameter  was  between  4.0  and  5.0,  and  for  M  —  8  and  16,  it  was  between 
5.0  and  6.0.  For  each  of  these  cases,  1000  widths  falling  between  the  stated  lower  and  upper 
bounds  were  used. 

11.2  Complex- Valued  Constellations 

The  exhaustive  searches  for  complex-valued  random  variables  were  not  as  fruitful  as  those 
for  real- valued  random  variables.  In  particular,  the  same  alphabet  size  (M  =  4)  and  bound 
on  symbol  magnitude  ( K  =  3)  parameters  as  used  in  the  real-valued  search  were  used  in  the 
search  for  complex- valued  variables.  The  constraint  (21)  was  omitted  in  light  of  the  results 
of  Appendix  A,  and  the  probabilities,  are  constrained  to  be  equal.  The  results  are  listed 
in  Table  8.  This  random  variable  corresponds  to  a  QPSK  or  4QAM  signal.  The  result  is 
significant  because  it  implies  that  the  symbol  alphabet  size  M  is  not  large  enough  to  allow 
the  simultaneous  satisfaction  of  the  various  constraints  and  small  cumulant  value.  However, 
searching  for  larger-alphabet  random  variables  is  prohibitively  computationally  costly  with 
available  computing  resources.  As  shown  in  the  next  section,  there  do  indeed  exist  larger 
alphabet  complex- valued  random  variables  with  small  fourth-order  cumulants. 

11.3  Low-Probability-of-Intercept  Signals 

In  this  section,  some  of  the  distributions  that  were  found  by  the  aforementioned  numerical 
techniques  are  used  to  design  digital  QAM  signals.  The  strengths  of  the  regenerated  sine 
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waves  (both  pure  and  impure)  for  orders  two,  three,  and  four  are  measured  and  tabulated, 
and  are  compared  to  the  strengths  of  BPSK,  16QAM,  and  duobinary  signals. 

The  complex-valued  random  variables  that  are  needed  to  create  the  complex-valued  PAM 
signals  are  constructed  from  the  minimum-cumulant  real-valued  random  variables  as  de¬ 
scribed  in  Section  9.  That  is,  let  b  and  c  be  independent  real-valued  random  variables  with 
some  minimum-cumulant  distribution.  The  complex-valued  random  variable  a  =  (b  +  ic)/y/ 2 
is  used  to  create  the  PAM  signal,  which  is  the  baseband  representation  of  the  radio-frequency 
digital  QAM  signal.  The  pulse  function  p(t)  is  that  for  duobinary  signaling  [3].  This  choice, 
as  noted  previously,  eliminates  the  second-order  cyclostationarity  that  is  associated  with  the 
pulse  rate.  Two  such  signals,  one  corresponding  to  row  three  of  Table  6  and  the  other  corre¬ 
sponding  to  row  one  of  Table  7,  were  simulated  and  their  second-,  third-,  and  fourth-order 
cyclic  moments  and  cumulants  for  a  =  0  and  1/T0  were  measured.  The  former  signal  is  de¬ 
noted  by  5\  and  the  latter  by  ST  The  results  are  shown  in  Table  10.  Similar  measurements 
for  three  common  digital  communication  signals  are  included  for  comparison.  The  signals 
were  generated  using  independent  bit  sequences,  and  the  total  data  length  is  32768  samples, 
which  corresponds  to  4096  symbols.  It  can  be  seen  from  the  table  that  the  strengths  of 
the  cumulant  sine  waves  for  the  signals  Si  and  S2  are  very  much  smaller  than  those  for  the 
other  three  signals,  indicating  that  the  two  signals  have  great  potential  to  avoid  detection 
by  SOCS  and  cumulant-based  interception  techniques.  Measurements  of  the  cyclic  moments 
were  also  made.  The  results  are  shown  in  Table  11.  The  results  in  the  table  indicate  that 
the  designed  signals  also  have  potential  to  avoid  detection  by  moment-based  interception 
techniques. 

Finally,  Table  9  provides  one  means  for  coding  a  bit  stream  into  the  symbols  corre¬ 
sponding  to  Si.  The  idea  is  to  code  each  incoming  pair  of  bits  (which  are  equiprobable  by 
assumption)  into  one  of  several  symbols  according  to  a  particular  probability  rule.  Thus, 
each  dibit  is  not  coded  into  a  unique  symbol  (except  for  the  dibit  00),  but  each  received 
symbol  can  be  decoded  into  only  one  dibit.  In  the  table,  the  third  column  contains  the 
probabilities  that  each  symbol  will  be  chosen  as  the  transmitted  symbol,  given  that  the  dibit 
in  the  first  column  is  presented  to  the  coder.  It  is  evident  that  the  price  paid  for  low  prob¬ 
ability  of  intercept  is  an  increased  number  of  states  in  the  modulation,  which  may  in  turn 
necessitate  more  power  in  the  tranmitted  signal  for  accurate  demodulation. 


12  Conclusions 

This  report  documents  a  study  of  the  possibility  of  designing  low-probability-of-intercept  sig¬ 
nals  by  minimizing  the  cumulants  of  the  complex-envelope  representation  of  digital  QAM  sig¬ 
nals.  In  particular,  digital  QAM  signals  at  radio  frequencies  can  be  represented  by  complex¬ 
valued  pulse-amplitude-modulated  signals  at  baseband,  and  the  cumulants  of  the  discrete 
random  variable  used  to  modulate  the  pulses  in  this  complex-envelope  signal  are  studied  for 
common  modulation  formats.  These  cumulants  are  directly  proportional  to  the  cumulants 
of  the  radio  frequency  signal.  Discrete  random  variable  distributions  are  sought  that  possess 
minimum  cumulants  and,  therefore,  that  minimize  the  cumulants  of  the  passband  signal. 

Several  numerical  methods  are  used  to  find  such  random  variable  distributions  (for  real¬ 
valued  variables),  complex- valued  partial-response  digital  QAM  signals  are  constructed  from 
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these  distributions,  and  their  cumulants  are  measured.  These  cumulants  are  shown  to  be 
much  less  than  those  for  ordinary  digital  QAM  signals,  and  are  less  than  those  for  commonly 
used  partial-response  signals.  In  addition,  the  moments  of  these  signals  are  measured  and 
are  also  seen  to  be  small.  Thus,  these  newly  designed  signals  have  potential  as  covert  signals. 

Searching  directly  for  complex-valued  random  variable  distributions  was  found  to  be  too 
computationally  costly  for  cases  of  interest.  In  particular,  the  size  of  the  symbol  constellation 
was  limited  to  four.  For  this  constellation  size,  the  minimum  fourth-order  cumulant  was 
found  to  be  equal  to  —1,  and  the  corresponding  constellation  is  that  for  QPSK  signaling.  This 
is  not  a  particularly  small  cumulant,  since  several  constellations  that  were  obtained  by  using 
real-valued  random  variables  yielded  fourth-order  cumulants  of  about  —0.01  and  —0.001,  but 
for  a  constellation  size  of  sixteen.  There  is  the  possibility  that  better  designs  exist,  but  they 
cannot  be  found  in  a  reasonable  amount  of  time  without  substantially  increased  computing 
speed. 

A  numerical  method  based  on  approximating  the  continuous  Gaussian  probability  density 
function  with  a  discrete  random  variable  was  found  to  be  computationally  simple,  but  yields 
random  variables  whose  cumulants  that  are  larger  than  those  found  by  exhaustive  search. 
If  distributions  with  larger  alphabets  than  sixteen  are  desired,  it  is  recommended  that  the 
exhaustive  search  method  be  used. 


PDF 

Size 

Constellation  Points  (Values  of  am) 

BPSK 

2 

±1 

4-level  PAM 

4 

(±1,  ±3)  x  75/5 

8-level  PAM 

8 

±1,  ±3,  ±5,  ±7  x  l/x/2l 

QPSK 

4 

±1,  ±z 

8PSK 

8 

±y/2/2  ±  y/2/2i,  ±1  ±  * 

16PSK 

16 

exp(z’27rg/16),  q  =  0,  •  •  • ,  15 

32PSK 

32 

exp(i27rq/32),q  =  0,  •  •  • ,  31 

64PSK 

64 

exp(z'27rg/64),  q  =  0,  •  •  • ,  63 

4QAM 

4 

~T^j2j2±vJW2 

8QAM 

n 

(±1,  ±z,  ±1  ±  i )  x  \/2/3 

16QAM 

16 

(±1  ±  i,  ±1  ±  *’/3,  ±1/3  ±  z/3,  ±1/3  ±  ?)  x  3/y/lO 

64QAM 

64 

(±<3,/7  ±  ir/7)  x  ^7/6,  for  q,  r  =  1, 3, 5,  7 

Table  1:  Numerical  values  of  the  points  in  the  symbol  constellations  that  correspond  to  some 
common  digital  communication  signals. 
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A  Cumulants  of  a  Constructed  Complex  Variable 

Suppose  b  and  c  are  independent  identically  distributed  real-valued  random  variables  such 

that  the  following  hold: 


E[b\ 

=  E[c]  =  0 

(24) 

E[b 2] 

=  E[S]  =  1 

(25) 

m 

0 

II 

”0 

II 

(26) 

E[b 4] 

=  E[c4]  =  C  +  3. 

(27) 

Find  the  cumulants  of  the  random  variable  a  =  b  +  ic. 


E[a]  =  CXfl  =  C*hl  =  E[b\  +  iE[c]  =  0 
£[|a|2]  =  62,1  =  E[b2  +  c2]  =  1  +  1  =  2 


15 


E[a2}  = 
E{a3}  = 
E[a2a }  = 

E[a4}  = 


E[a3a*}  = 

E[\a\4}  = 


^2,0  —  C^2 

E[b 2  -c2  +  2 ibc\  =  1  -  1  +  0  =  0 
C3fi  =  Cl3  =  E[(b  +  icf] 

E[b3  +  3i62c  +  36(ic)2  +  (ic)3]  =  0 
C^=Cl2  =  EHb+icf(b-ic)] 
£[i>3]  +  i£[c3]  +  i£[62c]  +  Elbe2]  =  0 
C4,0  =  C\A 
E[(b  +  ic)4] 


E 


bn(ic)4~n 


E[c4  +  b4-  6c2 b2]  =  2 C 

C4ti  =  C43  =  E[(b  +  ic)3(b  —  ic)} 

E[b4}  -  E[c4}  =  0 

C*4,2  +  2£7[|a  |2]2  =  C42  +  8 

E[b4  +  c4  +  262c2]  =  2  E[c4}  +  2 


C4, 2  =  2  C. 
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Constellation 

(PDF) 

M/C 

Order  n 

2 

4 

6 

8 

10 

BPSK 

M 

1.0 

1.0 

1.0 

1.0 

1.0 

4-level  PAM 

M 

1.0 

1.64 

2.92 

5.25 

9.45 

8-level  PAM 

M 

1.0 

1.76 

3.62 

7.92 

17.9 

BPSK 

C 

1.0 

-2.0 

16.0 

-272.0 

7936.0 

4-level  PAM 

c 

1.0 

-1.36 

8.32 

-111.85 

2603.14 

8-level  PAM 

c 

1.0 

-1.24 

7.19 

-92.0 

2039.5 

Table  2:  The  moments  (M)  and  cumulants  (C)  for  the  real-valued  signal  constellations 
shown  in  Table  1.  The  constellations  are  scaled  so  that  their  variance  is  one.  Moments  and 
cumulants  for  odd  values  of  n  are  equal  to  zero. 


PDF 

Order  of  Moment  ( n ) 

4 

6 

8  10 

Number  of  conjugated  variables  (m) 

KH 

2 

1,5 

3 

0,8 

ia 

4 

3,7 

5 

QPSK 

1.0 

1.0 

1.0 

1.0 

1.0 

| 

1.0 

1.0 

1.0 

1.0 

8PSK 

0.0 

1.0 

0.0 

DU 

1.0 

0.0 

mi 

MB 

0.0 

1.0 

16PSK 

0.0 

0.0 

1.0 

0.0 

0.0 

0.0 

■S3 

0.0 

BSM 

0.0 

0.0 

1.0 

0.0 

0.0 

1.0 

ua 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

1.0 

4QAM 

-1.0 

1.0 

-1.0 

1.0 

1.0 

-1.0 

1.0 

1.0 

8QAM 

-0.67 

1.68 

-1.48 

1.68 

16QAM 

-0.68 

-1.32 

1.96 

2.2 

-2.48 

3.12 

4.3 

64QAM 

-0.62 

-1.3 

2.23 

1.91 

-2.76 

3.96 

_ 

4.48 

-5.94 

7.58 

Table  3:  The  moments  for  the  complex-valued  signal  constellations  shown  in  Table  1,  All 
constellations  have  a  mean  of  zero  and  second-order  moment  (variance)  of  1.  The  moments 
for  values  of  m  not  shown  in  the  table  are  zero,  as  are  moments  for  odd  orders  n. 
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PDF 

Order  of  Cumulant  (n) 

4 

6 

8  10 

Number  of  conjugated  variables  (m) 

EO 

2 

ESI 

3 

0,8 

2,6 

4 

1,9 

3,7 

5 

QPSK 

1.0 

-1.0 

ESI 

-34.0 

34.0 

-34.0 

-496.0 

496.0 

8PSK 

0.0 

-1.0 

0.0 

4.0 

1.0 

0.0 

-33.0 

-8.0 

0.0 

456.0 

16PSK 

■mi 

-1.0 

0.0 

4.0 

0.0 

0.0 

-33.0 

0.0 

0.0 

456.0 

32PSK 

0.0 

-1.0 

ESI 

4.0 

0.0 

0.0 

-33.0 

0.0 

0.0 

456.0 

64PSK 

0.0 

-1.0 

0.0 

4.0 

0.0 

0.0 

-33.0 

WEM 

0.0 

456.0 

4QAM 

-1.0 

-1.0 

4.0 

4.0 

-34.0 

496 

496 

496 

8QAM 

-0.67 

-0.89 

2.30 

3.33 

-13.9 

-17.9 

-26.3 

178 

245 

352 

16QAM 

-0.68 

2.08 

-14.0 

-14.0 

162.7 

162.7 

64QAM 

1.8 

1.8 

-11.5 

-11.5 

-11.5 

127.5 

Table  4:  The  cumulants  for  the  complex- valued  signal  constellations  shown  in  Table  1.  All 
constellations  have  a  mean  of  zero  and  second-order  cumulant  (variance)  of  1.  The  cumulants 
for  values  of  m  not  shown  in  the  table  are  zero,  as  are  the  cumulants  for  odd  orders  n. 


n 

m 

Modulation  Type 

BPSK 

QPSK 

8PSK 

MPSK 
(M  >  16) 

APK 

k/T0±2fc 


k/T0 


k/T0 


k/T0  ±  4  / 


k/T0 


k/T0 


k/T0  ±  4  fc 


k/T0 


k/To  ±  4  /< 


k/T0 


k/To  ±  8 fc 


k/T0  ±  4/c  0 


k/To  k/To 


k/To  ±  8 fc  k/To  ±  8 fc 


10 

0,10 

k/To  ±  10  fc 

0 

0 

0 

10 

1,9 

k/To  ±  8/c 

k/To  ±  8/c 

k/T0  ±  8/c 

0 

10 

3,7 

k/To  ±  4/c 

k/T0  ±  4/c 

0 

0 

10 

k/To 

k/To 

k/To 

k/T0  ±  4  fc 


k/To 


k/To  ±  8 fc 


k/To  ±  Afc 


k/To 


k/To  ±  8/c 


k/To  ±  4  fc 


k/To 


Table  5:  Potential  cycle  frequencies  for  the  analytic  signals  corresponding  to  digital  QAM 
signals.  There  are  no  cycle  frequencies  for  the  values  of  m  that  are  not  shown  in  the  table, 
nor  for  odd  values  of  n.  The  +  sign  is  associated  with  the  first  value  of  m.  The  range  of  k 
for  a  given  signal  depends  on  the  excess  bandwidth. 
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5a,  5p 

Symbols  and  Probabilities 

Ca,4 

0.25,0.1 

2.0 

0.1 

-2.0 

0.1 

0.5 

0.4 

-0.5 

0.4 

0.25 

0.1, 0.1 

2.0 

0.1 

-2.0 

0.1 

0.5 

0.4 

-0.5 

0.4 

0.25 

0.05,0.1 

1.6 

0.2 

0.3 

0.2 

-0.35 

0.5 

-0.028 

0.05,0.05 

ieSi 

-0.85 

0.05 

-1.45 

0.25 

0.45 

0.65 

-0.00175 

Equiprobable  constraint 

0.05,  N/A  1.4 

0.25 

-1.4 

0.25 

0.15 

0.25 

-0.15 

0.25 

-1.03 

Table  6:  Results  of  a  numerical  search  for  real- valued  distributions  with  minimum  the  fourth- 
order  cumulant. 


M 

Symbols  and  Probabilities 

Ca,  4 

-1.75 

0.122 

1.75 

0.122 

-0.582 

0.378 

0.582 

0.378 

-0.65 

±2.29 

0.025 

±1.64 

0.070 

±0.98 

0.161 

±0.33 

0.244 

-0.32 

m 

±2.74 

0.005 

±2.37 

0.009 

±2.01 

0.020 

±1.64 

0.038 

-0.14 

IE3EH 

0.065 

±0.91 

0.096 

±0.55 

0.125 

±0.18 

0.143 

Table  7:  Results  of  approximating  the  Gaussian  distribution  to  obtain  real- valued  distribu¬ 
tions  with  minimum  the  fourth-order  cumulant. 


5a 

Symbols 

Ca,4 

0.5 

-1.0  ±  iO.O 

1.0  ±  iO.O 

0.0  +  i 

0.0  -  i 

-1.0 

0.25 

1.0  ±  iO.O 

0.0  ±  i 

0.0  -  i 

-1.0 

Table  8:  Results  of  a  numerical  search  for  complex-valued  distributions  with  minimum  the 
fourth-order  cumulant.  The  symbols  are  constrained  to  have  equal  probabilities. 
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Bit  Pair 

Constellation  Point 

Coding  Probability 

00 

(-0.35,  -0.35) 

1 

01 

(-2.0,  -2.0) 

0.04 

(1.6,  1.6) 

0.16 

(1.6,  -0.35) 

0.40 

(0.3,  -0.35) 

0.40 

10 

(-2.0,  1.6) 

0.08 

(1.6,  0.3) 

0.16 

(0.3,  1.6) 

0.16 

(-0.35,  -2.0) 

0.20 

(-0.35,  1.6) 

0.40 

11 

(-2.0,  0.3) 

0.08 

(-2.0,  -0.35) 

0.20 

(1.6,  -2.0) 

0.08 

(0.3,  -2.0) 

0.08 

(0.3,  0.3) 

0.16 

(-0.35,  0.3) 

0.40 

Table  9:  A  possible  coding  of  an  independent,  ecjuiprobable  bit  stream  into  the  symbols 
associated  with  signal  Si. 


n 

m 

a 

BPSK 

Duobinary 

16QAM 

5! 

S2 

2 

0 

2  fc 

0.0 

0.0 

-29.4 

-35.8 

-35.8 

2 

0 

2 fc  ±  1/To 

-9.7 

-120.0 

-37.4 

-120.0 

-120.0 

2 

1 

0 

0.0 

0.0 

0.0 

0.0 

0.0 

2 

1 

±1/T0 

-9.7 

-120.0 

-9.4 

-120.0 

-120.0 

3 

0 

3  fc 

-46.2 

-37.1 

-31.3 

-32.1 

-34.0 

3 

0 

3 fc  ±  1/To 

-56.6 

-62.1 

-39.2 

-59.5 

-55.2 

3 

1 

fc 

-46.2 

-37.1 

-35.6 

-40.9 

-30.3 

3 

1 

fc  ±  1/To 

-56.6 

-62.1 

-43.1 

-51.5 

-62.5 

4 

0 

4  fc 

6.0 

1.6 

-3.4 

-30.6 

-15.5 

4 

0 

4/c  ±  1/To 

-4.0 

-21.9 

-13.6 

-48.3 

-33.4 

4 

1 

2/e 

6.0 

1.6 

-27.4 

-25.5 

-33.5 

4 

1 

2/c  ±  1/To 

-4.0 

-21.9 

-38.6 

4 

2 

0 

6.0 

1.6 

-3.3 

-40.4 

4 

2 

±1/T0 

-4.0 

-21.9 

-12.2 

-46.8 

Table  10:  Peak  pure  (cumulant)  sine-wave  strengths  (in  dB)  for  various  complex-valued  dig¬ 
ital  communication  signals  for  orders  n  of  2,  3,  and  4,  and  for  various  numbers  of  conjugated 
factors  m.  S\  is  the  signal  obtained  by  using  a  minimum-cumulant  real-valued  distribution, 
and  S2  is  the  signal  obtained  by  using  a  minimum-cumulant  real-valued  distribution  obtained 
by  Gaussian  approximation. 
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a 

m 

a 

BPSI< 

Duobinary 

16QAM 

54 

S2 

B 

0 

2  fc 

0.0 

0.0 

-30.6 

-28.6 

-33.2 

B 

0 

2/c  ±  1/To 

-9.2 

-120.0 

-37.6 

-120.0 

-120.0 

B 

1 

0 

0.0 

0.0 

0.0 

0.0 

0.0 

B 

1 

±1/T0 

-9.2 

-120.0 

... 

-9.5 

-120.0 

-120.0 

B 

0 

3  fc 

-40.2 

-31.5 

-35.4 

-31.0 

-25.9 

B 

3 fc  ±  1/To 

iHS&Q 

-56.1 

-40.0 

-50.5 

-53.6 

B 

1 

fc 

-40.2 

-31.5 

-37.8 

-32.2 

-33.7 

B 

1 

fc  ±  1/To 

-49.8 

-56.1 

-43.0 

-53.8 

-51.1 

4 

0 

4  fe 

0.0 

5.3 

-3.3 

-32.3 

-15.1 

4 

0 

4 fc  ±  1/To 

-10.1 

-20.6 

-13.5 

-47.2 

-32.6 

4 

1 

2/c 

0.0 

5.3 

-36.4 

-17.4 

-26.1 

4 

1 

2/c  ±  1/To 

-10.1 

-20.6 

-45.4 

-43.6 

-40.5 

4 

2 

0 

0.0 

5.3 

2.5 

6.1 

5.2 

4 

2 

±1/T0 

-10.1 

-20.6 

-7.9 

-50.6 

-38.6 

Table  11:  Peak  impure  (moment)  sine-wave  strengths  (in  dB)  for  various  complex- valued  dig¬ 
ital  communication  signals  for  orders  n  of  2,  3,  and  4,  and  for  various  numbers  of  conjugated 
factors  m.  Si  is  the  signal  obtained  by  using  a  minimum-cumulant  real-valued  distribution, 
and  S2  is  the  signal  obtained  by  using  a  minimum-cumulant  real- valued  distribution  obtained 
by  Gaussian  approximation. 
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Signal  Detection  and  Sorting  Using  Cyclic  Cumulants 
in  the  General  Search  Algorithm 


by 


Chad  M.  Spooner  and  William  A.  Gardner 


Abstract 

The  problem  of  determining  the  number  of  cyclostationary  signals  that  are  present 
(if  any)  in  a  given  data  record  is  considered  in  this  report.  The  signals  can  be  com¬ 
pletely  spectrally  and  temporally  overlapping.  Nothing  is  assumed  about  the  signals 
except  that  they  exhibit  cyclostationarity  of  some  order.  The  detection  algorithm, 
called  the  general  search  algorithm ,  can  be  interpreted  as  a  blind  order- recursive  cyclic- 
cumulant  estimator,  and  is  tested  on  a  variety  of  signal  records  containing  one,  two, 
and  three  cochannel  signals.  The  results  indicate  that  the  algorithm  works  well  and 
merits  further  study  and  refinement. 
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1  Introduction 


When  confronted  with  the  problem  of  intercepting  signals  that  do  not  originate  from  friendly 
transmitters,  one  of  the  first  pieces  of  information  required  is  whether  or  not  a  signal  is 
present  in  the  given  data  record.  This  signal  detection  problem  becomes  a  bit  more  general — 
and  more  difficult — if  one  allows  for  the  possibility  of  the  presence  of  multiple  signals  in  the 
data.  In  this  report,  the  detection,  enumeration,  and  characterization  of  multiple  spectrally 
and  temporally  overlapping  unknown  signals  is  studied  from  the  point  of  view  of  exploiting 
the  structure  of  the  statistics  of  communication  signals.  In  particular,  many  (if  not  most) 
communication  signals  are  more  accurately  modeled  as  cyclostationary  signals  rather  than  as 
stationary  signals.  It  is  possible  to  exploit  the  properties  of  the  statistics  of  cyclostationary 
signals  to  successfully  perform  the  tasks  of  detection,  extraction,  direction-finding,  signal 
classification,  and  others  even  when  the  data  contains  multiple  temporally  and  spectrally 
overlapping  signals  [1]— [21] ,  [30,  31]. 

This  report  documents  an  initial  inquiry  into  the  possibility  of  using  the  higher-order 
moments  and  cumulants  of  a  cyclostationary  data  record  to  determine  the  number  of  signals 
present  in  the  record  and  to  estimate  some  of  their  modulation  parameters.  The  signals 
of  interest  herein  are  continuous-phase  frequency-shift  keying  (CPFSK),  phase-shift  keying 
(PSK)  and,  more  generally,  digital  quadrature-amplitude  modulation  (QAM),  although  the 
basic  search  algorithm  will  work  for  any  cyclostationary  signal.  It  is  the  processing  of  the 
output  of  the  search  algorithm  that  must  be  tailored  to  more  specific  signal  classes. 

A  typical  scenario  of  interest  is  described  in  the  following.  Suppose  a  given  data  record 
contains  a  binary  PSK  (BPSK)  signal,  a  quaternary  PSK  (QPSK)  signal,  and  an  offset  QPSK 
(OQPSK)  signal  as  well  as  additive  white  Gaussian  noise  (WGN).  This  complex-valued  data 
record  is  obtained  by  downconverting  and  sampling  the  original  radio-frequency  data.  The 
carrier  frequencies  of  the  signals  are  sufficiently  closely  spaced  so  as  to  preclude  linear  time- 
invariant  filtering  to  separate  the  signals.  The  processor  does  not  know  that  these  signals  are 
present  in  the  data  record.  The  goal  of  the  processing  is  to  determine  the  number  of  signals 
that  are  present  and  to  estimate  their  keying  rates  and,  possibly,  their  carrier  offsets  (the 
differences  between  their  carrier  frequencies  and  the  frequency  used  for  downconversion). 
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The  key  to  the  method  that  is  documented  in  this  report  is  that  the  signals  are  statistically 
independent,  which  implies  that  the  nth-order  temporal  cumulant  function  (TCF)  for  the 
data  is  asymptotically  equal  to  the  sum  of  the  nth-order  temporal  cumulant  functions  for  the 
individual  signals.  Each  of  these  cumulant  functions  is  a  periodic  or  polyperiodic  function  of 
time,  which  can  be  expressed  in  a  Fourier  series.  The  frequencies  of  the  Fourier  components 
are  related  to  the  modulation  parameters  (the  keying  rate  and  carrier  offset)  of  the  signals. 
If  the  Fourier  frequencies  for  each  individual  signal  are  distinct,  then  the  Fourier  components 
(sine  waves)  of  each  signal’s  ?rth-order  temporal  cumulant  function  can  be  estimated  from  the 
polyperiodic  nth-order  temporal  cumulant  function  of  the  data.  Thus,  using  the  cumulant 
provides  for  a  kind  of  signal  selectivity  in  the  processing.  Alternatively,  nth-order  temporal 
moment  functions  can  often  be  used  to  estimate  the  frequencies  of  the  Fourier  components 
of  TCFs,  but  these  parameters  are  not  signal  selective  in  general.  That  is,  the  amplitude 
and  phase  of  the  sine-wave  components  of  the  periodic  or  polyperiodic  nth-order  temporal 
moment  function  are  functions  of  all  signals  that  are  present  in  the  data,  but  only  some  of 
the  frequencies  are  functions  of  all  the  signals.  Thus,  temporal  moments  can  be  analyzed 
for  the  purpose  of  detection,  but  are  not  useful  for  classification  unless  there  is  only  one 
signal  present  (and  no  noise).  A  companion  report  presents  some  ideas  for  using  cumulants 
to  perform  signal-selective  modulation  classification. 

The  remainder  of  the  report  is  organized  as  follows.  In  Section  2  the  problem  of  interest 
is  stated  and  the  reasons  that  it  is  difficult  to  solve  are  discussed.  In  Section  3  the  proposed 
solution  is  described  in  detail,  and  some  illustrative  examples  are  provided  in  Section  4.  Sec¬ 
tion  5  presents  the  results  of  several  computer  simulations  for  signals  of  practical  interest, 
and  conclusions  are  drawn  in  Section  6.  The  theory  of  the  higher-order  statistics  of  cyclo¬ 
stationary  signals  that  forms  the  basis  for  the  approach  taken  here  is  reviewed  in  Appendix 
A. 

2  The  General  Search  Problem 

In  this  report,  we  are  interested  in  solving  the  general  search  problem.  The  general  search 
problem  is  defined  as  the  problem  of  detecting  the  presence  of  all  cyclostationary  signals 
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present  in  a  given  finite-length  data  record.  The  problem  is  difficult  in  that  there  is  little 
prior  knowledge  with  which  to  devise  an  algorithmic  solution.  Classical  decision  theory  is 
not  tractable  because  no  particular  probabilistic  models  for  the  signals  are  assumed  (ex¬ 
cept  that  they  are  cyclostationary),  the  signals  are  random  with  unknown  power  levels  and 
modulation  types,  and  the  noise  is  arbitrary.  If  the  signals  occupy  disjoint  portions  of  the 
spectrum,  then  energy  detection  and  filtering  could  be  used  to  detect  and  sort  the  signals, 
but  if  the  signals  are  spectrally  overlapping,  then  both  energy  detection  and  filtering  are 
ineffective.  The  kind  of  processing  necessary  is  that  for  which  the  output  of  the  processing 
can  be  effectively  partitioned  into  subgroups,  each  of  which  can  be  associated  with  one  and 
only  one  signal  in  the  data.  Since  linear  time-invariant  processing  is  not  effective  (filter¬ 
ing  cannot  separate  the  signals)  and  second-order  time-invariant  statistical  processing  is  also 
not  effective  (energy  detection),  both  nonlinear  and  time- varying  statistical  signal  processing 
are  possibilities.  As  it  turns  out,  cumulants  of  cyclostationary  signals  are  polyperiodically 
time- varying  parameters  obtained  by  nonlinear  operations. 

The  general  approach  is  motivated  by  the  two  assumptions  that  are  made:  (i)  the  signals 
are  cyclostationary  and  (ii)  the  signals  are  statistically  independent.  In  words,  and  at  the 
risk  of  oversimplifying,  the  approach  consists  of  estimating  the  amplitude,  frequency,  and 
phase  of  all  finite-strength  additive  sine-wave  components  that  appear  in  the  output  of 
specific  polynomial  transformations  of  the  data  (TCFs).  The  sine- wave  frequencies  are  then 
grouped  by  exploiting  harmonic  relationships. 

The  nth-order  TCF  can  be  analytically  calculated  for  a  large  class  of  communication 
signals.  This  class  is  defined  by  its  complex-envelope  representation,  which  is  a  complex¬ 
valued  pulse-amplitude-modulated  signal  with  independent,  identically  distributed  symbols 
and  arbitrary  pulse  function.  This  class  includes  as  special  cases  PSK,  00K,  and  M ary- 
QAM.  For  each  order  n  and  choice  of  the  number  of  optional  conjugations  m,  the  nth-order 
TCF  can  be  written  as  a  Fourier  series  in  which  each  sine-wave  component  has  complex¬ 
valued  amplitude  given  by  the  nth-order  cyclic  temporal  cumulant  function  (CTCF).  Thus, 
it  is  possible  to  list  the  amplitudes,  frequencies,  and  phases  for  each  order  for  signals  of 
this  class.  A  partial  list  of  potential  cycle  frequencies  is  given  in  Table  1.  As  evidenced 
by  the  table,  distinct  modulation  types  give  rise  to  distinct  patterns  (over  n  and  m)  of 
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cumulant  sine-wave  frequencies,  even  if  the  second-order  sine-wave  frequencies  are  identical. 
This  observation  can  be  used  to  devise  classification  methods  based  on  estimates  of  TCFs 
and  their  constituent  sine-wave  components.  For  the  purposes  of  detection  and  sorting,  we 
need  only  focus  on  the  case  in  which  the  number  of  conjugated  factors  m  is  half  the  order  n. 
A  consequence  of  this  restriction  is  that  only  cumulant  sine-waves  with  frequencies  equal  to 
harmonics  of  the  symbol  rates  of  the  signals  in  the  data  are  used,  and  these  frequencies  are 
enough  to  sort  and  detect  provided  that  the  symbol  rates  are  distinct.  If  the  symbol  rates  are 
not  distinct,  or  if  the  modulation  types  need  to  be  recognized,  then  other  choices  of  numbers 
of  conjugated  factors  must  also  be  included.  This  restriction  n  =  2m  also  substantially  eases 
the  computations  necessary  to  perform  the  detection  and  sorting  tasks  when  compared  to 
performing  these  tasks  based  on  the  combined  results  of  other  choices  of  m. 

The  algorithm  that  estimates  the  TCFs  (and,  thereby,  their  Fourier  components)  is  called 
the  general  search  algorithm  (GSA),  and  the  algorithm  that  groups  the  output  of  the  GSA 
is  called  the  grouping  algorithm  (GA).  These  algorithms  are  described  next. 

3  The  General  Search  and  Grouping  Algorithms 

The  approach  taken  to  solving  the  general  search  problem  consists  of  estimating  the  cycle 
frequencies  of  the  data  for  nonlinear  processing  of  various  orders  [5,  10,  8,  11].  In  order 
to  associate  the  resulting  cycle  frequency  estimates  with  specific  signals  in  the  data,  it  is 
required  to  estimate  cumulant  sine-wave  frequencies  rather  than  moment  cycle  frequencies. 
This  is  because  moment  cycle  frequencies  can  consist  of  sums  and  differences  of  the  cycle 
frequencies  for  various  distinct  signals  and  are,  therefore,  not  associated  with  any  particular 
signal  in  the  data,  and  because  the  strengths  of  the  desired  sine  waves  are  functions  of  all 
signals  in  the  data. 

Let  N  be  the  maximum  order  of  nonlinearity  that  is  to  be  used  for  processing.  The  goal 
of  the  processing  is  to  produce  a  list  of  cumulant  cycle  frequencies  {/?„}  for  each  value  of  n 
from  1  to  N.  The  list  {/?„}  characterizes  the  detectable  cyclostationarity  of  order  n  (and 
only  n )  that  is  associated  with  x(t )  because  it  is  not  contaminated  by  entries  that  are  due 
to  lower-order  sine  wave  interactions.  To  accomplish  this  task,  we  estimate  the  temporal 
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cumulant  function  (TCF)  for  x(t )  for  each  order  n.  From  this  estimate,  the  cycle  frequencies 
{/?„},  which  are  needed  for  the  estimate  of  the  TCF  for  order  n  +  1,  can  be  estimated. 

This  approach  is  justified  by  the  well-known  fact  that  the  periodogram  is  the  optimal 
estimator  of  the  frequency  of  a  sine  wave  in  white  Gaussian  noise.  Although  we  are  not 
primarily  interested  in  the  case  of  white  Gaussian  noise,  our  algorithm  essentially  estimates 
the  frequencies  of  a  set  of  relatively  strong  sine  waves  in  noise  and,  therefore,  the  periodogram 
is  of  interest.  In  addition,  the  algorithm  implements  an  estimator  of  the  pure-sine-waves 
function  for  each  order  n,  which  we  have  shown  is  the  signal-selective  function  of  interest  in 
applications.  More  explicitly,  the  general  search  problem  can  be  tackled  using  the  following 
general  search  algorithm  (GSA): 


0  Let  n  =  1,  fix  N  >  1,  denote  the  data  by  x(t),  0  <  t  <  T, 

choose  N  delays  T\,  •  •  •  Tyv,  and  choose  m  optional  conjugations. 

1  Compute  C'x{t,r)n  =  [n"=i  xH](t  +  Ti)]  “  E  p*  [llj=i  Cx{t,TUj)n]  for  r  =  [n 

2  Compute  Y(f)  =  FFT*  {<%(*,  r)n} 

3  Threshold  detect  the  bins  of  Y  to  find  {/?„} 

4  Compute  the  CTCFs  Cf"(r)n  =  (C,x(t,T)ne~i2^t)T 

5  Compute  the  TCF  Cx(t,r)n  =  C^{r)nei2^ 

6  n  — )■  n  +  1;  ifn^A^  then  go  to  1,  else  stop. 


V 


The  operation  of  the  GSA  and  the  notation  that  is  introduced  above  are  explained  in 
detail  next. 

In  Step  0,  the  maximum  order  of  nonlinearity  to  be  considered  is  fixed  at  N  >  1,  the 
N  delays  to  be  used  are  chosen,  the  m  optional  conjugations  are  chosen  for  each  processing 
order  n  <  N ,  and  the  processing  order  n  is  initialized  to  1. 

In  Step  1,  a  pre-estimate  of  the  nth-order  TCF  for  order  n  is  obtained  by  subtracting  from 
the  nth-order  delay  product  F[j=i  +  t3)  the  products  of  lower-order  TCFs  estimated 

in  previous  iterations  of  the  algorithm.  For  n  =  1,  there  are  no  previous  iterations,  so  the 
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first-order  pre-estimate  of  the  TCF  is  set  equal  to  the  first-order  lag  product  itself,  which  is 
just  the  data  x(t  +  Ti).  For  n  =  2,  the  product  of  the  first-order  TCF  estimates  for  each  of 
the  selected  lags  T\  and  r2  are  subtracted  from  the  second-order  lag  product.  This  removes 
from  consideration  any  sine  waves  in  the  second-order  lag  product  that  result  from  products 
of  first-order  sine  waves.  For  n  >  2,  the  sum  of  products  of  lower-order  TCFs  is  determined 
by  the  set  Pn,  which  is  the  set  of  distinct  partitions  of  the  set  of  indices  {1, 2,  •  •  • ,  n}.  This 
set  is  described  in  Appendix  A. 

In  Step  2,  the  pre-estimate  of  the  TCF  obtained  in  Step  1  is  Fourier  transformed  in  the 
t  variable  in  order  to  determine  its  sine-wave  components. 

In  Step  3,  the  values  of  this  transformed  TCF  pre-estimate  are  compared  to  a  threshold. 
The  locations  in  /  of  the  values  of  the  transformed  pre-estimate  that  exceed  the  threshold 
are  declared  to  be  cycle  frequencies  {/?„}. 

In  Step  4,  the  estimated  cycle  frequencies  are  used  to  compute  estimates  of  the  cyclic 
temporal  cumulant  functions  (CTCFs),  which  are  the  Fourier  coefficients  of  the  TCF  es¬ 
timates.  If  the  Fourier  transform  in  Step  2  has  length  equal  to  the  total  amount  of  data 
available  (T),  then  the  CTCFs  are  already  computed  in  Step  2,  and  do  not  need  to  be 
computed  again.  To  handle  the  case  in  which  the  cycle  frequencies  do  not  lie  on  bin-center 
frequencies,  interpolation  techniques  are  used  to  estimate  the  frequency  of  the  sine  wave, 
and  its  magnitude  and  phase  can  then  be  estimated  by  direct  computation  of  the  discrete 
Fourier  transform.  If  two  adjacent  bins  are  declared  to  correspond  to  sine-wave  frequencies, 
then  this  interpolation  is  done.  That  is,  if  two  adjacent  bins  have  large  magnitudes,  then  it  is 
assumed  that  a  single  sine  wave  with  frequency  somewhere  between  the  two  bin  frequencies 
is  responsible  for  this  energy. 

In  Step  5,  the  estimated  cycle  frequencies  and  CTCFs  are  combined  to  obtain  an  estimate 
of  the  TCF  that  replaces  the  pre-estimate  obtained  in  Step  1. 

Finally,  in  Step  6  the  order  of  processing  n  is  incremented  and  tested  against  its  maximum 
allowed  value  N.  If  n  is  less  than  or  equal  to  N  then  the  algorithm  returns  to  Step  1. 
Otherwise,  processing  is  terminated. 

Step  1  is  the  crucial  step.  Cycle  frequencies  could  be  estimated  by  Fourier  transforming 
the  lag  product  itself  and  thresholding  its  bins,  but  the  resulting  list  of  cycle  frequencies  would 
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contain  entries  due  to  interactions  among  the  distinct  signals  in  the  data.  By  subtracting 
the  particular  sum  of  products  of  lower-order  TCFs  in  step  1,  these  false  cycle  frequencies 
are  removed  from  consideration. 

The  output  of  the  GSA  is  a  sequence  of  lists  that  are  indexed  by  order.  Each  list  entry 
contains  two  elements.  The  first  is  the  cycle-frequency  estimate,  usually  denoted  by  a  or 
/?,  and  the  second  is  the  amplitude  of  the  sine  wave  with  frequency  a  for  the  appropriate 
order  (phase  information  is  suppressed  in  the  output  of  the  GSA,  but  retained  internally). 
Multiple  delay  sets  and  choices  of  conjugated  factors  can  be  accommodated  by  sequential 
runs  of  the  algorithm.  The  software  implementation  of  the  GSA  is  described  in  the  next 
section. 

3.1  The  GSA  Program 

The  GSA  program  takes  a  data  record  as  its  input  and  produces  cycle  frequency  and  cumulant 
estimates  for  specified  orders.  The  minimum  order  is  1  and  the  maximum  order  is  denoted 
by  A.  Typically,  the  minimum  order  is  set  to  2,  and — as  explained  below — only  the  even 
orders  between  2  and  N  are  used.  The  other  inputs  to  the  GSA  program  are  a  set  of  delay 
vectors  of  dimension  A,  and  a  set  of  conjugation  flags  of  dimension  N.  The  GSA  program 
then  estimates  the  Ath-order  temporal  cumulant  function  (TCF)  for  each  delay  vector  for 
each  conjugation  set.  For  instance,  the  program  could  be  used  to  compute  the  fourth-order 
cumulant  of  the  input  data,  for  the  two  delay  vectors 

12  5  6 
0  6  9  10 

and  the  two  conjugation-flag  vectors 

0  0  0  0 
0  0  10 

where  0  means  do  not  conjugate  and  1  means  conjugate.  This  would  result  in  four  fourth- 
order  temporal  cumulant  estimates.  Two  of  these  estimates  are  fourth-order  estimates  with 
0  conjugations,  and  the  remaining  two  are  fourth-order  estimates  with  the  third  variable 
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conjugated.  The  former  are  called  (4,  0)  estimates  (n  =  4  and  m  —  0),  and  the  latter 
are  (4,  1)  estimates.  Because  the  GSA  algorithm  is  recursive,  second-order  cumulants  are 
estimated  in  order  to  estimate  the  fourth-order  cumulants  (if  desired,  the  first-  and  third- 
order  estimates  can  be  made  as  well,  but  this  is  often  unnecessary  because  these  cumulants 
are  zero  for  almost  all  communication  signals  [except  those  with  pilot  tones]).  Since  TCFs 
are  polyperiodic  functions,  an  individual  TCF  estimate  can  be  represented  by  a  collection 
of  ordered  triplets  where  the  first  element  is  a  sine-wave  frequency,  the  second  is  a  sine-wave 
amplitude,  and  the  third  is  a  sine-wave  phase.  For  the  purpose  of  detecting  and  sorting, 
the  phase  is  not  needed.  Thus,  the  phase  is  retained  while  the  program  is  running,  but 
only  the  frequency  and  magnitude  parameters  are  actually  output.  An  individual  datum  in 
the  output  of  the  GSA  program  is  the  quadruple  (n,m,o,G),  where  n  is  the  order,  m  is 
the  number  of  conjugations,  a  is  the  cycle  frequency  estimate,  and  C  is  the  estimate  of  the 
magnitude  of  the  (n,  m)  CTCF  for  frequency  a.  These  quadruplets  are  indexed  by  the  delay 
vectors  r. 

The  output  of  the  GSA  program  contains  data  of  two  fundamental  sorts:  Xu,  which  is 
the  set  of  upper  data,  and  Xl,  which  is  the  set  of  lower  data.  These  are  defined  as 

n  =  2  m  Xl, 

n  ^  2m  =£•  x  6  Xu- 

Elements  in  Xl  have  cycle  frequency  estimates  a  that  are  not  (typically)  related  to  carrier 
frequencies,  whereas  elements  in  Xu  have  cycle  frequency  estimates  that  are  related  to  carrier 
frequencies. 

The  goal  of  the  grouping  algorithm,  and  its  implementation,  the  grouping  program ,  is  to 
group  this  multidimensional  data  into  sets  such  that  each  set  corresponds  to  one  and  only 
one  signal  in  the  original  data  record. 

3.2  The  Grouping  Algorithm 

The  output  of  the  GSA  program  is  a  set  of  lists  that  are  indexed  by  the  order  of  processing. 
Visual  inspection  of  these  lists  is  difficult  and  is  not  particularly  revealing.  Because  the 
cycle  frequencies  of  communication  signals  (especially  digital  QAM  signals)  are  harmonically 


9 


related  and  appear  at  multiple  orders  of  processing  n,  if  there  is  a  signal  present  in  the 
data,  then  there  will  be  a  set  of  a  estimates  that  are  harmonically  related.  The  main  idea, 
behind  the  grouping  algorithm  (GA)  is  to  extract  these  cycle  frequency  estimates  and  group 
them  together.  The  other  cycle  frequency-estimates  in  the  output  of  the  GSA  program  are 
discarded. 

The  word  “cluster”  in  the  following  description  of  the  GA  refers  to  a  standard  unsupervis- 
ed-learning  partitioning  algorithm  [32].  This  algorithm  finds  a  collection  of  subsets  of  a  given 
set  such  that  a  certain  cost  function  related  to  the  sample  mean  and  variance  of  each  subset 
is  minimized.  That  is,  at  termination  this  cost  would  increase  if  any  element  of  one  of  the 
sets  is  removed  and  then  added  to  any  of  the  other  sets. 

The  grouping  algorithm  consists  of  the  following  steps: 

1.  Read  GSA  data:  xj  =  (n,  m,  a,C)j  for  j  =  1,  •  •  • ,  M,  where  1  <  n  <  N  and  0  <  m  <  N 
for  each  j. 

2.  Separate  the  data  into  two  sets:  Xu,  which  is  the  set  of  upper  data  (n  ^  2m),  and  Xl, 
which  is  the  set  of  lower  data  (n  =  2m). 

3.  Cluster  the  set  Xl  into  three  sets  based  on  the  value  of  C.  That  is,  find  a  three-set 
partition  of  Xl  such  that  the  data  with  the  largest  C  are  in  one  set,  those  with  the 
smallest  C  are  in  another  called  Xw  (the  w  stands  for  “weak”),  and  the  rest  are  in  a 
third. 

4.  Find  the  union  of  the  two  sets  with  largest  C;  call  this  set  A"s  (the  s  stands  for  “strong”). 

5.  Cluster  Xs  based  on  the  value  of  a.  This  results  in  many  sets,  each  of  which  contains 
only  elements  with  a  that  are  “close  together.” 

6.  Search  the  set  Xw  for  any  elements  with  a  that  are  harmonically  related  (integer 
multiples  or  divisors)  to  the  mean  of  any  of  the  sets  obtained  in  the  previous  step.  If 
any  are  found,  add  them  to  an  appropriate  set.  Discard  the  remaining  elements  of  Xw. 

7.  Recluster  Xs. 
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8.  Group  the  sets  obtained  by  the  previous  step.  This  results  in  a  group  of  sets.  Each 
of  the  groups  is  associated  with  a  unique  fundamental  frequency.  The  members  of  the 
groups  are  sets  that  are  assumed  to  contain  cycle  frequencies  that  correspond  to  the 
various  harmonics  of  this  fundamental.  Compute  the  harmonic  numbers  for  each  of 
these  sets. 

9.  Store  the  vector  of  fundamentals  (one  fundamental  for  each  group)  for  later  use. 

10.  Cluster  the  set  Xu  into  three  sets  based  on  the  value  of  C.  That  is,  find  a  three-set 
partition  of  Xu  such  that  the  data  with  the  largest  C  are  in  one  set,  those  with  the 
smallest  C  are  in  another  called  Xw,  and  the  rest  are  in  a  third. 

11.  Find  the  union  of  the  two  sets  with  largest  C;  call  this  set  Xs. 

12.  Cluster  Xs  based  on  the  value  of  a. 

13.  Search  the  set  Xw  for  any  elements  with  a  that  are  separated  by  a  multiple  of  one 
of  the  stored  fundamentals  from  the  mean  of  any  of  the  sets  obtained  in  the  previous 
step.  This  must  be  done  only  for  data  that  have  matching  n  and  m  values.  If  any  such 
elements  are  found,  add  them  to  an  appropriate  set.  Discard  the  remaining  elements 
of  AG,. 

14.  Recluster  Xs. 

15.  Using  the  stored  fundamentals,  associate  the  sets  obtained  in  the  previous  step  with 
a  group  obtained  in  step  8.  Thus,  upper  sets  are  associated  with  lower  sets  through  a 
fundamental  frequency. 

16.  Estimate  the  carrier  offset  for  the  upper  elements  in  each  group. 

17.  Compute  the  harmonic  of  each  upper  cluster  in  each  group  by  using  the  fundamental 
and  the  estimated  offset. 

18.  Splinter  each  set  in  the  following  way.  Form  a  separate  set  for  each  of  the  distinct  (n, 
m)  pairs  that  appear  in  the  set.  At  the  end  of  this  procedure,  every  set  will  correspond 
to  a  single  (n,  m)  pair  and  to  a  single  harmonic. 
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19.  Output  the  matrix  of  detected  harmonics  for  each  group.  The  rows  of  these  matrices 
correspond  to  the  harmonic  number,  and  the  columns  correspond  to  each  (n,  m)  pair 
associated  with  that  group.  The  value  of  an  element  of  the  matrix  is  the  maximum 
value  of  the  C  parameters  of  all  the  data  points  x j  contained  in  the  set  that  corresponds 
to  the  appropriate  group  and  (n,  m)  pair.  This  matrix  shall  be  called  a  feature  matrix. 

4  Illustration  of  the  General  Search  Algorithm 

In  this  section,  some  simple  examples  of  the  operation  of  the  GSA  are  presented.  These 
examples  are  meant  to  illustrate  the  properties  of  cumulants,  estimates  of  which  are  used 
to  solve  the  general  search  problem.  More  realistic  examples  are  presented  and  discussed  in 
Section  5,  where  the  properties  of  the  cumulant  are  more  difficult  to  appreciate  directly. 

Since  the  ?ith-order  temporal  cumulant  function  is  also  the  pure-sine-waves  function  [11], 
the  nth-order  TCF  for  any  polyperiodic  signal  is  equal  to  the  signal  itself  for  n  =  1  and  is 
identically  zero  for  n  >  1.  This  is  because  there  can  be  no  pure  nth-order  sine  waves  for 
n  >  1  for  such  a  signal — all  sine-wave  components  in  the  higher-order  delay  products  must 
result  from  products  of  first-order  sine  waves.  A  simple  example  of  such  a  signal  is  a  single 
sine  wave.  Other  examples  include  the  sum  of  a  finite  number  of  sine  waves  and  any  periodic 
function,  such  as  a  square  wave  (which  is  equal  to  a  sum  of  an  infinite  number  of  sine  waves). 

Consider  a  sine  wave  with  frequency  1/11  and  amplitude  1.0.  Estimates  of  the  nth-order 
TCFs  for  orders  one  through  four  obtained  by  processing  a  sampled  version  of  such  a  sine 
wave  are  shown  in  Figure  1.  As  is  evident  from  the  figure,  the  first-order  TCF  is  equal  to 
the  sine-wave  itself,  and  the  other  TCFs  are  zero.  Next,  consider  the  sum  of  two  sine  waves 
with  frequencies  1/11  and  1/16  and  amplitudes  1.0  and  2.0,  respectively.  Estimates  of  the 
nth-order  TCFs  for  this  signal  are  shown  in  Figure  2  for  orders  one  through  four.  Again, 
only  the  first-order  TCF  is  nonzero.  A  third  sine  wave  with  frequency  1/7  is  added  and  the 
measurements  repeated.  The  results  are  shown  in  Figure  3.  Finally,  a  binary  full-duty-cycle 
rectangular-pulse  PAM  signal  with  pulse  rate  1/23  is  added  to  a  sine  wave  with  frequency 
1/11  and  the  measurement  procedure  is  repeated.  The  results  are  shown  in  Figure  4.  The 
second-  and  fourth-order  TCF  for  the  PAM  signal  for  the  chosen  delay  sets  should  be  square 
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waves  with  48%  duty  cycles.  As  can  be  seen  in  the  figure,  the  sine  wave  does  not  affect  the 
estimation  of  the  two  TCFs  for  the  PAM  signal,  which  do  not  appear  as  square  waves  with 
sharp  corners  due  to  truncation  of  the  Fourier-series  representation  of  the  TCF. 

These  examples  illustrate  that  nth-order  cumulants  do  indeed  compute  pure  nth-order 
sine  waves.  Thus,  the  GSA,  which  is  based  on  TCFs,  holds  promise  for  computing  the  correct 
cycle  frequencies  and  corresponding  Fourier  magnitudes  for  each  signal  in  a  given  data  set. 
The  following  section  presents  examples  of  the  operation  of  the  GSA  for  signals  of  interest. 


5  Computer  Simulations 

The  first  example  shows  the  capabilities  of  the  GSA  and  GA  program  for  the  case  of  16,384 
samples  of  a  binary  rectangular-pulse  PAM  signal  (To  =  23)  in  a  small  amount  of  WGN. 
The  temporal  cumulants  of  order  2,  4,  6,  and  8  are  estimated  for  six  delay  sets, 

00000000 
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00000007 
00000009 
0000000  12 

and  9  choices  of  conjugated  factors  (applied  to  each  of  the  six  delay  sets), 
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Thirty  cycle  frequencies  per  cumulant  are  estimated  and  used  in  the  computation  of  sub¬ 
sequent  cumulants;  only  fifteen  cycle  frequencies  per  cumulant  are  actually  output.  This 
is  done  for  the  sake  of  constraining  the  amount  of  computation:  the  run  time  of  the  GA 
program  is  a  function  of  the  total  number  of  cycle  frequencies  output  by  the  GSA  program. 
In  this  case,  and  all  subsequent  cases,  the  candidate  symbol  interval  lengths  are  constrained 
to  be  greater  than  2.0  and  less  than  25.0.  The  search  grid  has  fineness  0.1.  The  remaining 
parameter  to  explain  is  the  interval  parameter,  which  specifies  the  number  of  intervals  in 
which  to  segment  the  transform  of  each  pre-estimate  of  the  TCF.  The  maxima  of  each  in¬ 
terval  are  found  sequentially.  The  intervals  parameter  is  set  to  thirty.  The  feature  matrix 
(explained  subsequently)  for  the  output  of  the  grouping  program  is  shown  in  Figure  5. 

The  feature  matrix  consists  of  the  detected  harmonics  versus  the  order /conjugation  pairs, 
and  the  brightness  of  the  squares  indicates  the  strength  of  the  detected  harmonic.  For  all 
the  feature  matrices  shown  in  this  report,  the  first  three  columns  corresponds  to  order  two, 
the  next  five  correspond  to  order  four,  the  next  seven  to  order  six,  and  the  remaining  nine 
to  order  eight.  For  a  single  order,  the  columns  start  with  m  =  0  and  end  with  m  =  n.  In 
addition,  only  the  positive  harmonics  are  shown.  If  a  negative  harmonic  is  detected  and  the 
corresponding  positive  harmonic  is  not  detected,  then  the  positive  harmonic  is  replaced  by 
the  negative.  This  effectively  reduces  the  number  of  rows  in  the  plot  while  representing  the 
same  information.  Finally,  the  matrix  entries  that  correspond  to  a  =  0  for  each  order  are 
always  set  to  zero  (these  are  the  matrix  entries  for  n  —  2m  and  k  =  0).  The  cyclic  cumulants 
for  these  cycle  frequencies  are  not  signal  selective — all  signals  in  the  data  can  contribute  to 
these  cyclic  cumulants.  Thus,  they  are  not  useful  for  defining  a  feature  that  depends  only 
on  a  single  signal. 

In  the  second  example,  the  measurement  above  is  repeated  for  the  case  of  a  duobinary 
PAM  signal,  which  exhibits  no  second-order  cyclostationarity.  The  result  is  shown  in  Figure 
6. 

The  next  example  shows  the  capabilities  of  the  programs  for  the  case  of  16,384  samples  of 
a  complex- valued  BPSK  signal  with  symbol  rate  1/23  and  frequency  offset  0.004.  The  same 
cumulants  are  measured  with  the  same  parameters  as  in  the  previous  case.  The  feature 
matrix  for  this  case  is  shown  in  Figure  7.  The  grouping  program  correctly  estimates  the 
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symbol  rate  and  carrier  offset.  The  same  procedure  is  used  to  estimate  and  group  the  cycle 
frequencies  for  QPSK,  8PSK,  and  OQPSK  signals.  The  results  are  shown  in  Figures  8,  9, 
and  10. 

By  examining  Figures  5-10,  and  reviewing  the  relevant  entries  in  Table  1,  it  can  be  seen 
that  the  observed  patterns  are  correct,  although  there  may  be  a  few  harmonics  missing  here 
and  there.  In  particular,  the  feature  matrix  for  BPSK  should  be  “full,”  that  is,  all  harmonics 
are  exhibited  by  the  signal  for  all  even  orders  and,  indeed,  the  feature  matrix  for  BPSK  does 
show  harmonics  for  every  order/conjugation  pair  that  is  output  by  the  GSA  program.  For 
QPSK,  only  the  (n,  m)  pairs  such  that  n  —  2m  is  a  nonzero  multiple  of  4  should  exhibit 
carrier-related  cyclostationarity  (see  Table  1),  and  for  8PSK,  only  the  (8,  0)  and  (8,  8) 
pairs  should  exhibit  carrier-related  cyclostationarity.  The  case  of  offset  (staggered)  QPSK 
is  particularly  interesting  because  its  higher-order  statistics  cannot  be  analyzed  by  using 
the  previously  derived  cumulant  formulas  for  complex  PAM  because  its  complex  envelope 
is  a  PAM  signal  that  does  not  have  independent  symbols.  Nevertheless,  the  second-order 
cyclostationarity  for  this  signal  is  known  (see  [15],  page  447),  and  it  matches  with  the  first 
three  columns  of  the  Figure  10.  This  example  shows  that  the  GSA/GA  programs  can  be 
used  not  only  for  detection  and  sorting,  but  also  as  research  tools  that  facilitate  a  numerical 
evaluation  of  the  cyclostationarity  of  a  given  data  record. 

The  next  set  of  examples  shows  the  capability  of  the  algorithms  to  sort  two  equipower 
signals.  The  same  parameters  are  sued  as  in  the  previous  cases.  The  first  example  shows 
the  capability  of  the  programs  to  detect  and  sort  two  real-valued  rectangular-pulse  PAM 
signals.  One  of  the  signals  has  a  symbol  rate  of  1/15,  and  the  other  has  a  symbol  rate 
of  1/23.  The  signals  have  equal  power  levels.  The  feature  matrices  are  shown  in  Figures 
11-12.  Some  harmonics  are  missing  because  the  number  of  cycle  frequency  estimates  output 
per  delay  set  and  conjugation  set  pair  is  not  increased  with  respect  to  the  measurements 
that  were  performed  on  single  signals.  Also,  the  bottom  row  of  the  feature  matrix  for  12  is 
missing  because  the  two  PAM  signals  have  the  same  “carrier  offset”  of  zero.  Only  one  of 
the  signals  can  claim  these  cycle  frequency  estimates,  and  in  this  case  it  is  the  first  signal. 
The  ability  to  correctly  sort  the  cycle  frequency  estimates  for  the  case  of  signals  with  shared 
cycle  frequencies  is  to  be  added  at  a  later  date. 
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BPSK  signals  are  considered  next.  One  of  the  signals  has  carrier  offset  of  -0.005  and  a 
symbol  rate  of  1/15,  and  the  other  has  offset  0.004  and  a  symbol  rate  of  1/23.  The  feature 
matrices  are  shown  in  Figures  13-14.  This  test  is  repeated  for  two  QPSK,  8PSK,  and  OQPSK 
signals  that  have  the  same  symbol  rates  and  carrier  frequencies  as  the  two  BPSK  signals. 
The  results  are  showm  in  Figures  15-16,  17-18,  and  19-20,  respectively. 

The  next  set  of  examples  shows  the  capability  of  the  algorithms  to  sort  three  equipower 
signals.  For  these  cases,  only  the  lower  cycle  frequencies  are  used  (those  for  n  —  2m).  This 
is  done  both  to  demonstrate  that  this  is  a  viable  approach,  and  to  reduce  the  computational 
requirements  of  the  experiments.  The  cases  of  rectangular-pulse  PAM,  BPSK,  QPSK,  8PSK, 
and  OQPSK  are  considered.  The  results  are  shown  in  Figures  21-33. 

6  Conclusions 

This  report  documents  an  initial  inquiry  into  the  possibility  of  using  higher-order  cyclosta- 
tionarity,  and  in  particular  higher-order  cyclic  cumulants,  to  detect  and  sort  an  unknown 
number  of  cyclostationary  signals  in  a  given  data  record.  The  nth-order  cyclic  cumulants 
of  complex- valued  data  record  (analytic  signal  representation  of  a  radio-frequency  signal) 
for  orders  2  through  8  are  estimated  and  grouped  such  that  the  groups  each  correspond  to 
one  and  only  one  signal  in  the  data.  The  method  was  qualitatively  tested  for  several  signal 
environments  consisting  of  one,  two,  and  three  cochannel  signals  and  was  seen  to  perform 
well. 


7  Research  and  Development 

The  following  list  of  tasks  are  suggested  by  the  research  that  is  documented  in  this  report: 

1.  Detailed  computer  simulations  for  signal  environments  of  interest, 

2.  Mathematical  analysis  of  the  higher-order  cyclostationarity  of  communication  signals 
that  are  not  well  modeled  as  complex-valued  PAM  signals  at  baseband,  such  as  CPFSK, 
FSK,  and  analog  FM, 
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3.  Cataloging  of  the  higher-order  cyclic  features  of  signals  of  interest, 

4.  Mathematical  characterization  of  the  quality  of  the  GSA’s  cyclic  cumulant  estimates 
as  a  function  of  SNR,  SIR,  and  collect  time, 

5.  Development  of  generalizations  of  the  GA  to  handle  the  case  in  which  signals  share 
certain  cyclostationarity  properties  (e.g.,  two  or  more  signals  with  the  same  symbol 
rate  but  distinct  carrier  frequencies), 

6.  Study  alternate  techniques  for  estimating  sine-wave  frequencies  (rather  than  using  the 
discrete  Fourier  transform). 

7.  Study  alternate  feature-extraction  scheme  that  uses  the  cycle  frequencies  computed  by 
the  GSA  and  grouped  by  the  GA  to  do  precise  cumulant  measurements. 


n 

m 

Modulation  Type 

BPSK 

QPSK 

8PSK 

M  PSK 
(M  >  16) 

APK 

2 

0,2 

k/T0  ±  2 fc 

0 

0 

0 

m 

2 

1 

k/T0 

k/T0 

k/T0 

k/T0 

k/To 

4 

k/T0  ±  4/c 

k/T0  ±  4/c 

0 

k/T0  ±  4/c 

4 

2 

k/T0 

k/T0 

k/T0 

k/T0 

k/T0 

n 

0,6 

k/T0  ±  6  fc 

0 

0 

aa—gni 

6 

1,5 

k/T0  ±  4/c 

k/T0  ±  4/c 

0 

k/T0  ±  4/c 

6 

3 

k/T0 

k/T0 

k/T0 

k/T0 

k/To 

8 

0,8 

?r 

H- 

00 

k/T0  ±  8  fc 

k/T0  ±  8/c 

0 

\wsn/KC9iim\ 

8 

2,6 

k/T0  ±  4/c 

k/T0  ±  4/c 

0 

0 

k/T0  ±  4/c 

8 

4 

k/To 

k/T0 

k/T0 

k/T0 

k/To 

0,10 

0 

0 

El 

1,9 

k/T0  ±  8  fc 

k/T0  ±  8  fc 

k/T0  ±  8/c 

0 

k/T0  ±  8/c 

10 

3,7 

k/T0  ±  4/c 

k/T0  ±  4/c 

0 

0 

k/To  ±  4/c 

10 

5 

k/T0 

k/T0 

k/T0 

k/T0 

k/To 

Table  1:  Cycle  frequencies  for  the  analytic  signals  corresponding  to  PSK  and  other  digital 
QAM  signals.  APK  signals  include  non-PSK  digital  QAM  signals  used  in  practice,  such  as 
4QAM,  8QAM,  etc.  There  are  no  cycle  frequencies  for  the  values  of  m  that  are  not  shown 
in  the  table,  nor  for  odd  values  of  n.  The  +  sign  is  associated  with  the  first  value  of  m. 
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Figure  1:  The  temporal  cumulant  functions 
for  a  single  sine  wave.  The  delays  are  all  zero 
for  each  order. 


Figure  2:  The  temporal  cumulant  functions 
for  the  sum  of  two  sine  waves.  The  delays  are 
all  zero  for  each  order. 
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Figure  3:  The  temporal  cumulant  functions 
for  the  sum  of  three  sine  waves.  The  delays 
are  all  zero  for  each  order. 


t 

Figure  4:  The  temporal  cumulant  functions 
for  the  sum  of  a  sine  wave  and  a  full-duty- 
cycle  rectangular-pulse  binary  PAM  signal. 
The  delay  sets  are  [0],  [0  11],  and  [0  11  0  0], 
for  orders  one,  two,  and  four,  respectively. 
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Figure  5:  Measured  feature  matrix  for 
rectangular-pulse  PAM. 


Figure  7:  Measured  feature  matrix  for 
complex-valued  BPSK. 
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Figure  6:  Measured  feature  matrix  for  duobi-  Figure  8:  Measured  feature  matrix  for 
nary  (partial-response)  PAM.  complex-valued  QPSK. 
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o  5  io  is  20  0  Figure  12:  Measured  feature  matrix  for  one  of 

the  groups  output  by  the  GSA/GA  programs. 
Figure  10:  Measured  feature  matrix  for  The  input  is  the  sum  of  two  rectangular-puse 


complex- valued  offset  QPSK. 


PAM  signals.  The  program  correctly  identi¬ 
fied  the  symbol  rate  as  1/15,  and  carrier  offset 
as  -0.005. 
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Figure  13:  Measured  feature  matrix  for  one  of  Figure  15:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs,  the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  two  BPSK  signals.  The  input  is  the  sum  of  two  QPSK  signals. 
The  program  correctly  identified  the  symbol  The  program  correctly  identified  the  symbol 
rate  as  1/23,  and  carrier  offset  as  0.004.  rate  1/23,  and  carrier  offset  as  0.004. 
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Figure  14:  Measured  feature  matrix  for  one  of  Figure  16:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs,  the  groups  output  by  the  GS A / G A  programs. 
The  input  is  the  sum  of  two  BPSK  signals.  The  input  is  the  sum  of  two  QPSK  signals. 
The  program  correctly  identified  the  symbol  The  program  correctly  identified  the  symbol 
rate  as  1/15,  and  carrier  offset  as  -0.005.  rate  as  1/15,  and  carrier  offset  as  -0.005. 
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Figure  17:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  two  8PSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/23,  but  the  carrier  offset  as  0.032. 


Figure  19:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  two  OQPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/23,  and  carrier  offset  as  0.004. 


Figure  18:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  two  8PSI<  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/15,  but  the  carrier  offset  as  0.059. 


Figure  20:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  two  OQPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/15,  and  carrier  offset  as  -0.005. 
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Figure  21:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  rectangular- 
pulse  PAM  signals.  The  program  correctly 
identified  the  symbol  rate  as  1/19. 


Figure  22:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  rectangular- 
puse  PAM  signals.  The  program  correctly 
identified  the  symbol  rate  as  1/15. 


Figure  23:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  rectangular- 
puse  PAM  signals.  The  program  correctly 
identified  the  symbol  rate  as  1/23. 


Figure  24:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  BPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/23. 
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Figure  25:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  BPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/15. 


Figure  27:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  QPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/23. 
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Figure  26:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  BPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/19. 


Figure  28:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  QPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/19. 
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Figure  29:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  QPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/15. 
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Figure  30:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  8PSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/23. 
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Figure  31:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  8PSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/19. 


Figure  32:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  8PSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/15. 
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Figure  33:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  OQPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/15. 
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Figure  35:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  OQPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/11.5. 


Figure  34:  Measured  feature  matrix  for  one  of 
the  groups  output  by  the  GSA/GA  programs. 
The  input  is  the  sum  of  three  OQPSK  signals. 
The  program  correctly  identified  the  symbol 
rate  as  1/9.5. 
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A  The  Theory  of  Higher-Order  Cyclostationarity 

In  this  appendix  we  define  the  temporal  and  spectral  moment  and  cumulant  functions  that 
form  the  basis  of  the  theory  of  IIOCS.  Then  we  give  a  brief  tutorial  explanation  of  how  cum- 
ulants  arise  as  the  solution  to  the  problem  of  generating  pure  nth-order  sine  waves.  Finally, 
we  explain  the  signal-selectivity  property  that  is  unique  to  the  cyclic  temporal  cumulants 
and  their  Fourier  transforms,  the  cyclic  polyspectra,  and  we  illustrate  the  parameters  using 
the  example  of  digital  QAM  signals. 

For  a  time-series  x(t)  for  — oo  <  t  <  o o,  we  define  the  nth-order  lag-product  time-series 
by 

n 

Lx(t,r)n=  ITX(*  +  U)>  C1) 

3=1 

where  r  =  and  [-]f  denotes  matrix  transposition.  The  cyclic  temporal  moment 

function  (CTMF)  of  order  n  is  defined  by  the  limiting  time  average 

K  (r)n  =  lim  —  /  Lr(t,T)ne-'2’°,dt  =  (L,(t,T)„e-2'’°,\,  (2) 

T-+  OO  1  J-T/2  ' 

which  is  simply  the  Fourier  coefficient  associated  with  the  component  el2nat  in  the  time- 
series  Lx{t,r)n.  It  can  be  seen  that  the  CTMF  is  a  Fourier  coefficient  of  a  moment  function 
because  the  nth-order  fraction-of-time  probabilistic  moment  (the  temporal  moment  function 
[TMF])  associated  with  the  lag  product  Lx(t,r)n  is  given  by 

Rx(t,T)n=  E^{Lx(t,r)n}  (3) 

where  E^  {•}  =  ((•)e-t27raQ  el2nat,  and  can  be  expressed  as  [15,  13] 

Rx{t,r)n=  £  Rax(r)ne'2™\  (4) 

oefa) 
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where  the  sum  is  over  all  real  numbers  a,  called  nth-order  cycle  frequencies ,  for  which 
RaAT)n±  o.  In  (3),£<“>{  •}  is  the  temporal  expectation  operation  (or  the  sine-wave  extrac¬ 
tion  operation).  The  functions  (2)  and  (3)  exist  and  are  well-behaved  for  appropriate  mod¬ 
els  of  many  time-series  including  amplitude  modulated  (AM),  pulse-amplitude-modulated 
(PAM),  phase-shift-keyed  modulated  (PSK),  and  digital  quadrature  AM  (QAM)  signals 
[10,  4],  and  others. 

The  temporal  cumulant  function  (TCF)  of  order  n  for  the  set  of  time-series  translates 
{ x(t  +  Tj )}”_j  is  defined  by 


Cx(t,r)n  =  E 

p=M  L, 


p 

k{p)  U  i 


j=  1 


(5) 


which  is  completely  analogous  to  its  stochastic-process  counterpart  in  [22].  The  sum  in  (5)  is 
over  all  distinct  partitions  Pn  of  the  set  of  indices  {1, 2,  •  •  • ,  n},  where  each  partition  {vif}Vk=\ 
has  p  elements,  1  <  p  <  n,  and  k(p)  =  (— l)p~1(p  —  1)!.  The  vector  rUj  is  the  vector  of  n3 
lags  with  indices  in  the  set  Uj. 

The  cyclic  temporal  cumulant  function  (CTCF)  of  order  n  is  the  Fourier  coefficient  of  the 
TCF: 

C£(t)„£  (6) 

The  set  of  real  numbers  {0}  for  which  Cf(r)n  ^  0  is  called  the  set  of  pure  nth-order  cycle 
frequencies ,  for  reasons  that  will  become  clear  subsequently.  Combining  (3)-(6)  reveals  that 
the  CTCF  is  given  by  the  following  explicit  function  of  lower-order  CTMFs: 


Cx  (r)n  =  E 
P 


k(p )  E  (  II 

a 1 1 =0  v=i 


(7) 


where  1  =  [1  •  •  •  1]*,  and  a  =  [ai  ■  •  •  ap]h  The  CTCF  was  originally  derived  in  [2]  as  the 
solution  to  the  problem  of  removing  from  the  Fourier  coefficient  R^(t)„  all  contributions 
from  Fourier  coefficients  i?rJ(r^)nj  of  lower  order.  This  is  equivalent  to  removing  from  the 
finite-strength  additive  sine-wave  component  of  frequency  0  in  the  lag  product  time-series 
Lx(t:  r)„  all  contributions  from  products  of  sine-wave  components  in  lag  products  Lx(t ,  TVj)n3 
of  lower  order — whose  frequencies  sum  to  0 — that  can  be  obtained  by  factoring  Lx(t,r)n. 
By  using  the  relationship  between  moments  and  cumulants[27],  the  CTMF  can  be  expressed 
in  terms  of  CTCFs  of  order  n  and  lower: 


K(r)n  =  E 

p 


e  n  Gx3  (r^ ) nj 

L/3fl=aJ'=1 


(8) 


The  CTMF  and  the  CTCF  are  not  in  general  integrable  due  to  the  presence  of  sinusoidal 
components  in  r.  These  components  formally  result  in  Dirac  deltas  in  the  n-dimensional  Fou¬ 
rier  transform  of  the  CTCF.  However,  a  reduced-dimension  version  of  the  CTCF  is  absolutely 
integrable  for  many  time-series  of  interest  and,  therefore,  it  is  strictly  Fourier  transformable 
[2].  The  reduced-dimension  CTCF  (RD-CTCF)  is  simply  the  CTCF  associated  with  the  n 
variables  {a ft  +  Tj)}”_ x  with  rn  =  0.  The  RD-CTCF  is  denoted  by 


Cf(u)n  =  Cf(r)n  for  =  ii(  and  t„  =  0, 


(9) 
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where  u  is  (n  —  l)-dimensional  vector  [rii  •  •  •  The  (n  —  l)-dimensional  Fourier  transform 

of  (9)  is  denoted  by  Pf(f')n: 


ps(f)n  =  r  ■■■  r  c?(u)n  e-a*u'f’du,  (io) 

J— oo  J—oo 

where  /'  =  [/i  •  •  •  /»_ i]f. 

The  spectral  moment,  spectral  cumulant,  and  cyclic  polyspectrum  are  defined  as  follows. 
Consider  the  n  complex-demodulate  time-series  Xr(t,fj )  for  j  =  1  , •••,n,  associated  with 
narrow  bandpass  filtered  versions  of  x(t ),  where 

rt+T/2 

XT(t,f)  =  /  x(v)e~'2*Jv  dv.  (11) 

Jt-T/2 

The  limit  as  T  — >  oo  of  the  limiting  time-average  of  the  product  of  these  spectral  components 
is  called  the  spectral  moment  function  (SMF)  of  order  n 

«/)»=  Km  (nA'T(i, /,-)),  (12) 

and  it  can  be  shown  that  Dirac  deltas  can  be  factored  out  as  follows: 

&(/)»  =  ££?(/V(/'i-«).  (i3) 

a 

where  S(-)  is  the  Dirac  delta  function.  However,  the  factor  S°(f)n  contains  additional  Dirac 
deltas  for  many  signals  and  n  >  2  (e.g.,  BPSK  and  n  =  4). 

The  spectral  cumulant  function  (SCF)  of  order  n  is  given  by 


Px(f)n 


£ 


hp)  n 


j=i 


(14) 


where  /  is  the  vector  of  n3  frequencies  with  subscripts  in  the  set  i/j,  and  it  follows  from 
(13)  that  Dirac  deltas  can  again  be  factored  out: 


W)»  =  £Pf(/')n<5(/'l-/3).  (15) 

0 


Analogous  to  the  definition  of  the  cyclic  spectrum  (and  power  spectrum)  in  [15],  the 
factor  Pf(f')n  is  defined  to  be  the  cyclic  polyspectrum  (CP)  and  is  given  explicitly  by 


P,6(f)n  =  £ 


P-1 


Ctfl=/3 


3=1 


(16) 


As  first  shown  in  [2],  the  CP  is  the  (n  —  l)-dimensional  Fourier  transform  of  the  RD-CTCF 
Cf(u)n  (cf.  (10)).  This  is  a  generalization  of  the  Wiener  relation  between  the  power  spec¬ 
trum  and  autocorrelation  from  second-order  stationary  time-series  (cf.  [15])  to  nth-order 
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cyclostationary  time-series.  Within  the  stochastic-process  framework  of  generally  nonsta¬ 
tionary  processes,  it  should  be  called  the  cyclic  Shiryaev-Kolmogorov  relation  [22],  which  is 
the  generalization  of  the  Wiener-Khinchin  relation  (cf.  [15]).  Because  the  RD-CTCF  can  be 
absolutely  integrable,  the  CP — unlike  the  SMF — contains  no  Dirac  deltas.  That  is,  all  Dirac 
deltas  present  in  the  individual  terms  in  (16)  cancel. 

Let  us  now  see  how  the  cumulant  arises  as  the  solution  to  the  problem  of  generating 
pure  nth-order  sine  waves.  For  low  orders  n,  it  is  easy  to  mathematically  characterize  a 
pure  nth-order  sine  wave  in  a  way  that  matches  our  intuition.  For  notational  simplicity, 
we  choose  the  case  of  no  conjugated  factors  in  (1).  For  n  =  1,  the  moment  sine  waves 
(the  sine-wave  components  of  the  polyperiodic  TMF)  are,  by  definition,  pure  lst-order  sine 
waves.  For  n  =  2,  all  products  of  lst-order  moment  sine  waves  can  be  subtracted  from  the 
2nd-order  moment  sine  waves  to  obtain  the  pure  2nd-order  sine  waves,  which  are  denoted 
by  <xr(*,T1,T2) 2  , 

M*,ti,t2)2  =  E{a]{x(t  +  n)x(t  +  t2)}  -  E{a}  {x{t  +  n)}  E{a}  {x(t  +  t2)} 

=  Rx(t,r)2  -  Rx(t,n)i  Rx(t, t2)i. 

There  are  several  interesting  points  to  be  made  concerning  pure  2nd-order  sine  waves: 
(i)  since  Rx(t,  r,)i,  i  =  1,2,  and  Rx{t,r)2  are  1st-  and  2nd-order  moments  respectively,  then 
ax(t,  Ti,  r2)2  is  a  temporal  covariance  function;  (ii)  if  Rx(t ,  t)x  =  0,  then  there  are  no  lower- 
order  sine  waves,  and  the  2nd-order  moment  sine  waves  are  equal  to  the  pure  2nd-order 
sine  waves;  (in)  if  the  variables  x(t  +  rx)  and  x(t  +  r2)  are  statistically  independent  (in  the 
temporal  sense  [15,  13]),  then  E ^  {x(t  -t-  Tx)r(i  +  r2)}  —  {x(t  +  Ti)}  E {x(t  +  r2)} 
and  therefore  crx(t,  rx,  r2)2  =  0,  that  is,  there  is  no  pure  2nd-order  sine  wave  for  this  particular 
pair  of  delays  rx  and  r2.  This  latter  point  shows  that  the  set  of  pure  cycle  frequencies  {/?} 
can  be  considerably  different  from  the  set  of  impure  (moment)  cycle  frequencies  {a}. 

The  pure  3rd-order  sine  waves  are  obtained  next.  From  the  3rd-order  moment  sine  waves, 
we  want  to  subtract  each  possible  product  of  lower-order  sine  waves,  but  only  once  each. 
Thus,  we  subtract  products  of  pure  2nd-order  and  pure  lst-order  sine  waves  from  the  3rd- 
order  moment  sine  waves,  rather  than  subtracting  products  of  1st-  and  2nd-order  moment 
sine  waves: 


crx{t,  t)3  =  E{a}  I JJ  x(t  +  Tj)  |  -  ax(t,  Tj,  t2)2  crx(t,  r3)  1  -  ax(t,  n,  r3)2  ax(t,  r2)i 

Crx{t,T2,T3)2  0’x(t,T\'ji  C7x{t^T 1)1  O x ( t ,  T2 ) X  CrI(t,T3)x  . 


Observe  that  there  are  no  other  possible  products  of  pure  lower-order  sine  waves.  The  terms 
in  the  sum  of  products  that  are  subtracted  can  be  enumerated  by  considering  the  distinct 
partitions  of  the  index  set  {1,2,3}.  A  partition  of  a  set  G  is  a  collection  of  p  subsets  of  G, 
with  the  following  properties:  G  =  Uj=i  vj  and  Vj  0^  =  0  for  j  ^  k.  The  set  P3 
of  distinct  partitions  of  {1,2,3}  is  given  by 

P  =  1:  {1,2,3} 

P  =  2:  {1,2},  {3}  {1,3},  {2}  {2,3},  {1} 

P  =  3:  {!},  {2},  {3}. 
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Thus,  we  can  express  the  pure  3rd-order  sine  wave  ax(t7r) 3  as  a  sum  over  the  elements  of 

P3: 

p 

ax{t,T)3  =  Rx(t,r) 3-J2  I]>*(^K  ,  (17) 

P3  b=1 

p/i 

where  t„.  is  a  subset  of  {Tj}l=1  with  elements  having  indices  in  vj,  and  nj  is  the  number  of 
elements  in  uj.  Notice  that,  as  in  the  case  of  n  =  2,  if  the  lst-order  moment  sine  waves  are 
zero,  then  the  3rd-order  moment  sine  waves  are  equal  to  the  pure  3rd-order  sine  waves.  In 
this  case,  there  are  no  products  of  lower-order  sine  waves  to  subtract  from  the  moment  sine 
waves. 

The  formula  for  the  pure  nth-order  sine  waves  is 


Tc(C  ^”)n  —  Rx{t,T)n  '/’  ] 


Pn 
Pit- 1 


JJ  Tr(C  )Uj 

j= 1 


(18) 


where  Pn  is  the  set  of  distinct  partitions  of  the  index  set  {1,2,  •  •  •  ,n}.  The  pure-sine-waves 
formula  (18)  gives  all  the  pure  nth-order  sine  waves  associated  with  the  delay  set  r.  A 
single  pure  nth-order  sine  wave  with  frequency  (3  can  be  selected  by  computing  the  Fourier 
coefficient: 

<xf  (r)n  ei2*0t  =  (crx(u,  r)n  e~^^)  ,  (19) 

and  can  be  expressed  in  terms  of  pure  lower-order  sine  waves  by  using  the  Fourier  series  for 
each  o-x(t,T„.)nj  in  (18): 


ax(t7w)k  =J2crxk(.w)kel27r0k\ 

pk 


W  =  [tUl  •  •  -tU/;]1, 


(20) 


where  the  sum  is  over  all  cycle  frequencies  j3k  of  order  k.  Thus,  the  strength  of  a  single  pure 
nth-order  sine  wave  is  given  by 


<?x(T)n  =  Rx(r)n  -  E 


Pn 

Pltl 


£ 

1 


n 

j= 1 


A 


(TiJj  )n] 


(21) 


where  (3  is  the  p-dimensional  vector  of  pure  cycle  frequencies  \(3\  ■  ■  ■  /3pY  and  1  is  the 
p-dimensional  vector  of  ones.  Hence,  the  pure-sine-wave  strength  <r^(r)n  is  given  by  the 
CTMF  Rx(r)n  with  all  products  of  pure  lower-order  sine-wave  strengths,  for  sine  waves 
whose  frequencies  sum  to  (3,  subtracted  out.  From  (8),  it  is  clear  that  the  pure-sine-wave 
strength  <rf(r)n  is  identical  to  the  CTCF  Cf(r)„. 

Finally,  let  us  see  how  the  CTCF  (and  the  CP)  is  signal  selective.  Let  the  signal  x(t) 
consist  of  the  sum  of  M  mutually  independent  signals  {ym(0}m=n 

M 

x(t)  =  H  (22) 

m= 1 
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Then  the  addition  rule  for  cumulants  [4,  23]  can  be  used  to  show  that  the  nth-order  TCF 
for  x(t)  is  the  sum  of  TCFs  for  {ym(t)}, 

M 

C,(t,T)»=ZC,m(t,T)„.  (23) 

771=1 

Thus,  the  pure  nth-order  sine  waves  in  the  delay  products  of  each  ym(t)  add  to  form  the 
pure  nth-order  sine  wave  in  the  delay  product  of  x(t): 

M 

C&T).  =  £  Cl(T)„.  (24) 

771=1 

The  TMF  does  not  in  general  obey  this  useful  cumulative  relation  (24).  That  is,  the  nth-order 
TMF  for  x{t)  is  not  the  sum  of  nth-order  TMFs  for  each  ym(t):  Rx(t,  r)n  /  Ylm= i  Rym(ti  r)n 
.  An  exception  is  the  case  of  zero-mean  signals  and  n  <  3,  for  which  moments  and  cumulants 
are  equal,  which  is  a  commonly  encountered  case  in  HOS. 

To  illustrate  how  (24)  can  be  applied  in  practice,  consider  the  situation  in  which  {ym(t)} 
represent  M  interfering  signals  that  overlap  in  time  and  frequency,  but  each  ym(t)  possesses 
some  distinct  nth-order  cycle  frequency,  say  j3m.  Then  it  follows  from  (24)  that 

Cf“M«  =  C*T(t)„,  (25) 

This  indicates  that  the  presence  or  absence  of  each  of  the  signals  ym(t)  can  be  detected 
by  measuring  (estimating)  the  CTCFs  of  x(t)  for  the  cycle  frequencies  {/3m},  and  that 
parameters  of  each  of  the  signals,  on  which  these  CTCFs  depend,  can  be  estimated.  As 
illustrated  in  [15,  12,  14]  for  second  order  and  in  [4]  for  higher  order,  this  signal-selectivity 
property  can  be  exploited  in  numerous  ways  to  accomplish  noise-and-interference-tolerant 
signal  detection  and  estimation. 

As  another  application,  let  M  =  2,  yi(t)  be  non-Gaussian,  and  j/2 (^)  be  Gaussian.  Then 
Cy2{t,r)n  =  0  for  n  >  3  and,  from  (23),  we  have  Cx(t,r)n  —  Cyi(t,T)n ,  n  >  3,  which 
indicates  the  detectability  of  yi(t)  with  no  knowledge  about  j/2(^)  except  that  it  is  Gaussian. 

As  an  example  of  the  cyclic  polyspectrum,  we  consider  the  class  of  digital  QAM  signals 
described  by  their  complex  envelopes,  which  can  be  expressed  as  x(t)  =  J2m=-oo  amP(t  — 
mTo  —  to),  where  p(t)  is  the  complex  pulse  and  {am}  is  the  complex  digital  data.  Assuming 
that  { am }  is  an  independent  sequence,  we  have  shown  [10,  5,  11]  that  the  CP  for  the  set  of 
variables  {od*b(£  -j-  Tj)}"=1  (with  rn  —  0),  where  (*)j  denotes  an  optional  conjugation  of  the 
jth  variable,  is  given  by 

%(f) n  =  -  l’/'])w"  n  D  =  k/To,  (26) 

10  j=  1 

where  (— )j  is  the  optional  minus  sign  corresponding  to  the  optional  conjugation  (*)j,  Ca,n 
is  the  cumulant  of  am,  and 

/OO 

p(t)e~W  it.  (27) 

-oo 

Thus,  the  CP  is  simply  a  scaled  product  of  n  pulse  transforms. 
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Abstract 

In  this  report  the  problem  of  determining  the  number  and  modulation  types  of 
cyclostationary  signals  that  are  present  (if  any)  in  a  given  data,  record  is  studied.  The 
signals  are  assumed  to  be  cochannel;  that  is,  the  signals  are  assumed  to  substantially 
overlap  in  time  and  frequency.  The  general  search  algorithm  is  used  to  blindly  estimate, 
in  an  order-recursive  procedure,  higher-order  cyclic  cumulants,  then  the  grouping  al¬ 
gorithm  is  used  to  associate  the  cyclic  cumulants  with  each  other  and,  therefore,  with 
signals  in  the  data,  and  finally  the  classification  algorithm  is  used  to  determine  the 
modulation  type  for  each  signal  based  on  the  grouped  cyclic  cumulant  estimates.  The 
GSA/GA  tandem  is  described  in  detail  and  applied  to  simulated  data  to  produce  a 
catalog  of  features  for  classification.  These  two  algorithms  are  applied  to  data  records 
that  contain  two  cochannel  signals  to  illustrate  the  fact  that  the  measurements  of  the 
features  are  signal  selective.  Finally,  the  CA  is  described  and  it  is  shown  that  cer¬ 
tain  signal  types  cannot  be  uniquely  classified  (given  the  cochannel  signal  assumption) 
without  using  cyclic  cumulants  of  order  larger  than  two. 
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1  Introduction 


There  are  several  situations  in  which  a  signal  reception  system  receives  multiple  signals  with 
unknown  modulation  types.  One  example  is  the  general  nuisance  of  cochannel  interference. 
In  this  situation,  the  signal  of  interest  is  received  along  with  other  signals  that  are  not  of 
interest,  and  about  which  nothing  is  known.  Another  example  is  when  the  radio  spectrum  is 
monitored  for  the  purpose  of  regulating  radio  transmissions  in  a  given  geographical  area.  In 
the  former  situation,  the  task  is  to  remove  the  effects  of  the  interfering  signals  on  the  signal 
of  interest,  and  in  the  latter  situation,  the  task  is  to  determine  the  types  of  signals  that 
are  received  and,  possibly,  to  demodulate  the  signals  to  determine  their  message  content.  A 
third  example  is  that  of  signal  search  for  interception,  in  which  the  goal  is  to  determine  as 
much  information  as  possible  about  all  the  unknown  signals  that  are  present  (including  the 
modulation  format,  location  and  type  of  emitter,  and,  possibly,  the  message). 

When  a  single  signal  is  received  in  a  small  amount  of  noise,  various  techniques  can 
be  used  to  make  a  decision  about  the  signal’s  modulation  type  [25]— [37].  For  example,  a 
phase  histogram  can  be  computed,  the  bandwidth  and  center  frequency  can  be  accurately 
estimated,  the  modulus  can  be  measured,  pilot  tones  can  be  detected,  etc.  Various  automatic 
modulation  classifiers  can  be  constructed  based  upon  these  and  other  measurements.  On  the 
other  hand,  when  multiple  signals  are  received  and  these  signals  are  substantially  spectrally 
and  temporally  overlapping,  or  when  the  noise  is  strong  with  respect  to  the  signal,  these 
techniques  tend  to  fail.  The  primary  reason  for  this  is  that  the  above  measurements  are 
not  signal  selective;  each  signal  (including  noise)  contributes  to  the  measurements,  thereby 
degrading  them,  and  reducing  the  efficacy  of  the  subsequent  classification  algorithm. 

In  this  report,  the  possibility  of  using  the  cyclostationarity  property  of  communication 
signals  to  perform  modulation  recognition  is  investigated.  In  particular,  the  second-  and 
higher-order  cyclostationarity  [14,  11,  4,  5]  of  the  signals  is  used  to  create  signal-selective 
features  for  modulation  classification  and  recognition.  Because  the  relevant  statistical  pa¬ 
rameters  associated  with  cyclostationarity  are  asymptotically  independent  of  stationary  noise 
and  interference,  so  are  the  derived  feature  vectors.  In  addition,  if  the  modulation  param¬ 
eters  of  the  received  signals  are  distinct  (e.g.,  distinct  baud  rates  and  carrier  frequencies), 
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then  the  feature  vector  for  a  particular  received  signal  will  be  asymptotically  independent  of 
all  other  received  signals. 

2  The  Classification  Problem 

Let  x(t)  denote  the  received  data,  which  consists  of  the  sum  of  M  statistically  independent 
signals  Sj(t)  together  with  additive  noise: 

M 

x{t)  =  ^2sj(t)  +  w(t)  0  <t<T.  (1) 

j= i 

M  is  assumed  to  be  unknown  and  w(t)  is  assumed  to  be  stationary  and  Gaussian.  The  signal 
x(t)  is  complex- valued,  and  it  is  assumed  that  it  is  obtained  by  in-phase  and  quadrature  sam¬ 
pling  of  the  downconverted  radio-frequency  signal.  Thus,  if  the  downconversion  frequency 
is  equal  to  the  carrier  frequency  of  the  signal  Sk(t),  then  this  signal  has  a  carrier  offset  of 
zero.  The  carrier  offsets  of  the  M  signals  are  assumed  to  be  sufficiently  closely  spaced  so  as 
to  preclude  the  use  of  filtering  to  separate  the  signals. 

The  classification  problem  considered  herein  is  to  determine  the  modulation  type  of  each 
of  the  signals  Sj(t).  The  first  part  of  this  problem  consists  of  determining  the  number  of 
cyclostationary  signals  that  are  present  in  the  data  record.  This  is  called  the  general  search 
problem.  Higher-order  cyclostationarity  is  used  to  devise  a  solution  to  the  general  search 
problem.  This  solution  provides  estimates  of  a  feature  vector  for  each  signal  in  the  data. 
Classification  can  then  be  accomplished  by  using  the  measured  feature  for  each  signal. 

The  main  idea  behind  the  use  of  higher-order  cyclostationarity  to  perform  signal  classifi¬ 
cation  is  that  of  sine-wave  generation  by  nonlinear  transformation.  When  a  cyclostationary 
signal  is  subjected  to  a  series  of  nonlinear  transformations,  a  series  of  sets  of  finite-strength 
additive  sine  waves  are  produced.  Distinct  modulation  formats  produce  distinct  sets  of 
sine  waves  if  the  series  of  nonlinear  transformations  is  chosen  correctly.  When  the  series  of 
transformations  in  applied  to  data  that  contains  multiple  signals  with  distinct  modulation 
parameters,  the  resulting  sine  waves  can  be  grouped  into  disjoint  subsets  such  that  each 
subset’s  strengths,  frequencies,  and  phases  depend  (asymptotically)  on  only  a  single  signal 
that  is  present  in  the  data.  The  particular  series  of  nonlinear  transformations  of  the  data 
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that  is  proposed  here  is  the  series  of  nth-order  temporal  cumulant  functions. 

The  remainder  of  this  report  is  organized  as  follows.  In  Section  3,  an  overview  of  the 
classification  procedure  is  presented.  In  Section  4,  the  general  search  algorithm  (GSA),  which 
is  used  to  compute  the  sequence  of  nth-order  temporal  cumulant  functions,  and  the  group¬ 
ing  algorithm  (GA),  which  is  used  to  associate  (group)  the  detected  sine-wave  frequencies 
obtained  by  the  GSA,  are  presented  and  explained.  In  Section  5,  the  feature  vectors  for 
classification  are  defined  in  terms  of  the  output  of  the  GA,  and  extensive  computer  simu¬ 
lation  examples  are  presented.  In  Section  6,  the  proposed  classification  algorithm  (CA)  is 
presented.  The  fundamentals  of  the  theory  of  the  higher-order  statistics  of  cyclostationary 
signals  are  provided  in  Appendix  A. 


3  Overview  of  Solution 

A  crude  outline  of  the  classification  procedure  is  given  in  the  following: 

•  Collect  data. 

•  Fix  maximum  order  N  of  nonlinearity  for  processing. 

•  Estimate  temporal  cumulant  functions  (TCFs)  for  orders  1  (or  2)  through  N  using  the 
general  search  algorithm. 

•  Extract  from  the  estimated  TCFs  a  set  of  TCFs  for  individual  signals  in  the  data  using 
the  grouping  algorithm. 

•  Compare  each  extracted  set  of  TCFs  to  known  sets  of  TCFs  for  signals  of  interest. 
Classify  the  modulation  type  of  the  signal  based  on  the  similarity  to  the  known  TCFs 
using  the  classification  algorithm. 

4  The  General  Search  and  Grouping  Algorithms 

The  approach  taken  to  solving  the  general  search  problem  consists  of  estimating  the  cycle 

frequencies  of  the  data  for  nonlinear  processing  of  various  orders  [5,  10,  8,  11].  In  order 
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to  associate  the  resulting  cycle  frequency  estimates  with  specific  signals  in  the  data,  it  is 
required  to  estimate  cumulant  sine-wave  frequencies  rather  than  moment  cycle  frequencies. 
This  is  because  moment  cycle  frequencies  can  consist  of  sums  and  differences  of  the  cycle 
frequencies  for  various  distinct  signals  and  are,  therefore,  not  associated  with  any  particular 
signal  in  the  data,  and  because  the  strengths  of  the  desired  sine  waves  are  functions  of  all 
signals  in  the  data. 

Let  N  be  the  maximum  order  of  nonlinearity  that  is  to  be  used  for  processing.  The  goal 
of  the  processing  is  to  produce  a  list  of  cumulant  cycle  frequencies  {/?„}  for  each  value  of  n 
from  1  to  N.  The  list  {/?„}  characterizes  the  detectable  cyclostationarity  of  order  n  (and 
only  n)  that  is  associated  with  x(t)  because  it  is  not  contaminated  by  entries  that  are  due 
to  lower-order  sine  wave  interactions.  To  accomplish  this  task,  we  estimate  the  temporal 
cumulant  function  (TCF)  for  x(t)  for  each  order  n.  From  this  estimate,  the  cycle  frequencies 
{/?„},  which  are  needed  for  the  estimate  of  the  TCF  for  order  n  +  1,  can  be  estimated. 

This  approach  is  justified  by  the  well-known  fact  that  the  periodogram  is  the  optimal 
estimator  of  the  frequency  of  a  sine  wave  in  white  Gaussian  noise.  Although  we  are  not 
primarily  interested  in  the  case  of  white  Gaussian  noise,  our  algorithm  essentially  estimates 
the  frequencies  of  a  set  of  relatively  strong  sine  waves  in  noise  and,  therefore,  the  periodogram 
is  of  interest.  In  addition,  the  algorithm  implements  an  estimator  of  the  pure-sine- waves 
function  for  each  order  n,  which  we  have  shown  is  the  signal-selective  function  of  interest  in 
applications.  More  explicitly,  the  general  search  problem  can  be  tackled  using  the  following 
general  search  algorithm  (GSA): 
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Let  n  —  1,  fix  N  >  1,  denote  the  data  by  x(t),  0  <  t  <  T, 
choose  N  delays  r1?  •  •  -tjv,  and  choose  m  optional  conjugations. 

Compute  C'x(t,  r)n  =  [n”=i  x^(t  +  r,-)]  -  £  [llj=i  Cx(t,  rUj)n]  for  r  =  [tj 

p#  i 

Compute  Y(f)  =  FFT<  {C'x(t,r)n} 

Threshold  detect  the  bins  of  Y  to  find  {(3n} 

Compute  the  CTCFs  Cf»(r)n  =  (£'(*,  r)n 
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5  Compute  the  TCF  Cx{t,r)n  =  E/?„  C’f"(r)ne,'2^"< 

6  n  — >  n  +  1;  if  n  <  N  then  go  to  1,  else  stop. 

The  operation  of  the  GSA  and  the  notation  that  is  introduced  above  are  explained  in 
detail  next. 

In  Step  0,  the  maximum  order  of  nonlinearity  to  be  considered  is  fixed  at  N  >  1,  the 
N  delays  to  be  used  are  chosen,  the  m  optional  conjugations  are  chosen  for  each  processing 
order  n  <  N,  and  the  processing  order  n  is  initialized  to  1. 

In  Step  1,  a  pre-estimate  of  the  nth-order  TCF  for  order  n  is  obtained  by  subtracting  from 
the  nth-order  delay  product  flLi  x^(t  +  t3)  the  products  of  lower-order  TCFs  estimated 
in  previous  iterations  of  the  algorithm.  For  n  =  1,  there  are  no  previous  iterations,  so  the 
first-order  pre-estimate  of  the  TCF  is  set  equal  to  the  first-order  lag  product  itself,  which  is 
just  the  data  x{t  +  iy).  For  n  =  2,  the  product  of  the  first-order  TCF  estimates  for  each  of 
the  selected  lags  Ti  and  r2  are  subtracted  from  the  second-order  lag  product.  This  removes 
from  consideration  any  sine  waves  in  the  second-order  lag  product  that  result  from  products 
of  first-order  sine  waves.  For  n  >  2,  the  sum  of  products  of  lower-order  TCFs  is  determined 
by  the  set  Pn,  which  is  the  set  of  distinct  partitions  of  the  set  of  indices  {1,2,  •  •  •  ,n}.  This 
set  is  described  in  Appendix  A. 

In  Step  2,  the  pre-estimate  of  the  TCF  obtained  in  Step  1  is  Fourier  transformed  in  the 
t  variable  in  order  to  determine  its  sine-wave  components. 

In  Step  3,  the  values  of  this  transformed  TCF  pre-estimate  are  compared  to  a  threshold. 
The  locations  in  /  of  the  values  of  the  transformed  pre-estimate  that  exceed  the  threshold 
are  declared  to  be  cycle  frequencies  {/ 3n }. 

In  Step  4,  the  estimated  cycle  frequencies  are  used  to  compute  estimates  of  the  cyclic 
temporal  cumulant  functions  (CTCFs),  which  are  the  Fourier  coefficients  of  the  TCF  es¬ 
timates.  If  the  Fourier  transform  in  Step  2  has  length  equal  to  the  total  amount  of  data 
available  (T),  then  the  CTCFs  are  already  computed  in  Step  2,  and  do  not  need  to  be 
computed  again.  To  handle  the  case  in  which  the  cycle  frequencies  do  not  lie  on  bin-center 
frequencies,  interpolation  techniques  are  used  to  estimate  the  frequency  of  the  sine  wave, 
and  its  magnitude  and  phase  can  then  be  estimated  by  direct  computation  of  the  discrete 
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Fourier  transform.  If  two  adjacent  bins  are  declared  to  correspond  to  sine-wave  frequencies, 
then  this  interpolation  is  done.  That  is,  if  two  adjacent  bins  have  large  magnitudes,  then  it  is 
assumed  that  a  single  sine  wave  with  frequency  somewhere  between  the  two  bin  frequencies 
is  responsible  for  this  energy. 

In  Step  5,  the  estimated  cycle  frequencies  and  CTCFs  are  combined  to  obtain  an  estimate 
of  the  TCF  that  replaces  the  pre-estimate  obtained  in  Step  1. 

Finally,  in  Step  6  the  order  of  processing  n  is  incremented  and  tested  against  its  maximum 
allowed  value  N.  If  n  is  less  than  or  equal  to  N  then  the  algorithm  returns  to  Step  1. 
Otherwise,  processing  is  terminated. 

Step  1  is  the  crucial  step.  Cycle  frequencies  could  be  estimated  by  Fourier  transforming 
the  lag  product  itself  and  thresholding  its  bins,  but  the  resulting  list  of  cycle  frequencies  would 
contain  entries  due  to  interactions  among  the  distinct  signals  in  the  data.  By  subtracting 
the  particular  sum  of  products  of  lower-order  TCFs  in  step  1,  these  false  cycle  frequencies 
are  removed  from  consideration. 

The  output  of  the  GSA  is  a  sequence  of  lists  that  are  indexed  by  order.  Each  list  entry 
contains  two  elements.  The  first  is  the  cycle-frequency  estimate,  usually  denoted  by  a  or 
/?,  and  the  second  is  the  amplitude  of  the  sine  wave  with  frequency  a  for  the  appropriate 
order  (phase  information  is  suppressed  in  the  output  of  the  GSA,  but  retained  internally). 
Multiple  delay  sets  and  choices  of  conjugated  factors  can  be  accommodated  by  sequential 
runs  of  the  algorithm.  The  software  implementation  of  the  GSA  is  described  in  the  next 
section. 

4.1  The  GSA  Program 

The  GSA  program  takes  a  data  record  as  its  input  and  produces  cycle  frequency  and  cumulant 
estimates  for  specified  orders.  The  minimum  order  is  1  and  the  maximum  order  is  denoted 
by  N.  Typically,  the  minimum  order  is  set  to  2,  and — as  explained  below — only  the  even 
orders  between  2  and  N  are  used.  The  other  inputs  to  the  GSA  program  are  a  set  of  delay 
vectors  of  dimension  Ar,  and  a  set  of  conjugation  flags  of  dimension  N.  The  GSA  program 
then  estimates  the  ATh-order  temporal  cumulant  function  (TCF)  for  each  delay  vector  for 
each  conjugation  set.  For  instance,  the  program  could  be  used  to  compute  the  fourth-order 
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cumulant  of  the  input  data  for  the  two  delay  vectors 


12  5  6 
0  6  9  10 

and  the  two  conjugation-flag  vectors 

0  0  0  0 
0  0  10 

where  0  means  do  not  conjugate  and  1  means  conjugate.  This  would  result  in  four  fourth- 
order  temporal  cumulant  estimates.  Two  of  these  estimates  are  fourth-order  estimates  with 
0  conjugations,  and  the  remaining  two  are  fourth-order  estimates  with  the  third  variable 
conjugated.  The  former  are  called  (4,  0)  estimates  (n  =  4  and  m  =  0),  and  the  latter 
are  (4,  1)  estimates.  Because  the  GSA  algorithm  is  recursive,  second-order  cumulants  are 
estimated  in  order  to  estimate  the  fourth-order  cumulants  (if  desired,  the  first-  and  third- 
order  estimates  can  be  made  as  well,  but  this  is  often  unnecessary  because  these  cumulants 
are  zero  for  almost  all  communication  signals  [except  those  with  pilot  tones]).  Since  TCFs 
are  polyperiodic  functions,  an  individual  TCF  estimate  can  be  represented  by  a  collection 
of  ordered  triplets  where  the  first  element  is  a  sine- wave  frequency,  the  second  is  a  sine-wave 
amplitude,  and  the  third  is  a  sine-wave  phase.  For  the  purpose  of  detecting  and  sorting, 
the  phase  is  not  needed.  Thus,  the  phase  is  retained  while  the  program  is  running,  but 
only  the  frequency  and  magnitude  parameters  are  actually  output.  An  individual  datum  in 
the  output  of  the  GSA  program  is  the  quadruple  (n,m,a,  C),  where  n  is  the  order,  m  is 
the  number  of  conjugations,  a  is  the  cycle  frequency  estimate,  and  C  is  the  estimate  of  the 
magnitude  of  the  (n,  m)  CTCF  for  frequency  a.  These  quadruplets  are  indexed  by  the  delay 
vectors  r. 

The  output  of  the  GSA  program  contains  data  of  two  fundamental  sorts:  Xu,  which  is 
the  set  of  upper  data,  and  Xl,  which  is  the  set  of  lower  data.  These  are  defined  as 

n  =  2 m  =>•  x  €  Xl, 

n  7^  2 m  =>  x  G  Xu- 
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Elements  in  Xl  have  cycle  frequency  estimates  a  that  are  not  (typically)  related  to  carrier 
frequencies,  whereas  elements  in  Xv  have  cycle  frequency  estimates  that  are  related  to  carrier 
frequencies. 

The  goal  of  the  grouping  algorithm ,  and  its  implementation,  the  grouping  program ,  is  to 
group  this  multidimensional  data  into  sets  such  that  each  set  corresponds  to  one  and  only 
one  signal  in  the  original  data  record. 

4.2  The  Grouping  Algorithm 

The  output  of  the  GSA  program  is  a  set  of  lists  that  are  indexed  by  the  order  of  processing. 
Visual  inspection  of  these  lists  is  difficult  and  is  not  particularly  revealing.  Because  the 
cycle  frequencies  of  communication  signals  (especially  digital  QAM  signals)  are  harmonically 
related  and  appear  at  multiple  orders  of  processing  n,  if  there  is  a  signal  present  in  the 
data,  then  there  will  be  a  set  of  a  estimates  that  are  harmonically  related.  The  main  idea 
behind  the  grouping  algorithm  (GA)  is  to  extract  these  cycle  frequency  estimates  and  group 
them  together.  The  other  cycle  frequency-estimates  in  the  output  of  the  GSA  program  are 
discarded. 

The  word  “cluster”  in  the  following  description  of  the  GA  refers  to  a  standard  unsupervis- 
ed-learning  partitioning  algorithm  [24],  This  algorithm  finds  a  collection  of  subsets  of  a  given 
set  such  that  a  certain  cost  function  related  to  the  sample  mean  and  variance  of  each  subset 
is  minimized.  That  is,  at  termination  this  cost  would  increase  if  any  element  of  one  of  the 
sets  is  removed  and  then  added  to  any  of  the  other  sets. 

The  grouping  algorithm  consists  of  the  following  steps: 

1.  Read  GSA  data:  Xj  =  ( n ,  m,  oc,C)j  for  j  =  1,  •  •  • ,  M,  where  1  <  n  <  N  and  0  <  m  <  N 
for  each  j. 

2.  Separate  the  data  into  two  sets:  Xu,  which  is  the  set  of  upper  data  (n  ^  2m),  and  XL , 
which  is  the  set  of  lower  data  (n  =  2m). 

3.  Cluster  the  set  Xl  into  three  sets  based  on  the  value  of  C.  That  is,  find  a  three-set 
partition  of  Xl  such  that  the  data  with  the  largest  C  are  in  one  set,  those  with  the 
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smallest  C  are  in  another  called  Xw  (the  w  stands  for  “weak”),  and  the  rest  are  in  a 

third. 

4.  Find  the  union  of  the  two  sets  with  largest  C;  call  this  set  Xs  (the  s  stands  for  “strong”). 

5.  Cluster  Xs  based  on  the  value  of  a.  This  results  in  many  sets,  each  of  which  contains 
only  elements  with  a  that  are  “close  together.” 

6.  Search  the  set  Xw  for  any  elements  with  a  that  are  harmonically  related  (integer 
multiples  or  divisors)  to  the  mean  of  any  of  the  sets  obtained  in  the  previous  step.  If 
any  are  found,  add  them  to  an  appropriate  set.  Discard  the  remaining  elements  of  Xw. 

7.  Recluster  Xs. 

8.  Group  the  sets  obtained  by  the  previous  step.  This  results  in  a  group  of  sets.  Each 
of  the  groups  is  associated  with  a  unique  fundamental  frequency.  The  members  of  the 
groups  are  sets  that  are  assumed  to  contain  cycle  frequencies  that  correspond  to  the 
various  harmonics  of  this  fundamental.  Compute  the  harmonic  numbers  for  each  of 
these  sets. 

9.  Store  the  vector  of  fundamentals  (one  fundamental  for  each  group)  for  later  use. 

10.  Cluster  the  set  Xu  into  three  sets  based  on  the  value  of  C.  That  is,  find  a  three-set 
partition  of  Xu  such  that  the  data  with  the  largest  C  are  in  one  set,  those  with  the 
smallest  C  are  in  another  called  Xw,  and  the  rest  are  in  a  third. 

11.  Find  the  union  of  the  two  sets  with  largest  C;  call  this  set  Xs. 

12.  Cluster  Xs  based  on  the  value  of  a. 

13.  Search  the  set  Xw  for  any  elements  with  a  that  are  separated  by  a  multiple  of  one 
of  the  stored  fundamentals  from  the  mean  of  any  of  the  sets  obtained  in  the  previous 
step.  This  must  be  done  only  for  data  that  have  matching  n  and  m  values.  If  any  such 
elements  are  found,  add  them  to  an  appropriate  set.  Discard  the  remaining  elements 
of  Xw. 
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14.  Recluster  Xs. 


15.  Using  the  stored  fundamentals,  associate  the  sets  obtained  in  the  previous  step  with 
a  group  obtained  in  step  8.  Thus,  upper  sets  are  associated  with  lower  sets  through  a 
fundamental  frequency. 

16.  Estimate  the  carrier  offset  for  the  upper  elements  in  each  group. 

17.  Compute  the  harmonic  of  each  upper  cluster  in  each  group  by  using  the  fundamental 
and  the  estimated  offset. 

18.  Splinter  each  set  in  the  following  way.  Form  a  separate  set  for  each  of  the  distinct  (n, 
m)  pairs  that  appear  in  the  set.  At  the  end  of  this  procedure,  every  set  will  correspond 
to  a  single  (n,  m)  pair  and  to  a  single  harmonic. 

19.  Output  the  matrix  of  detected  harmonics  for  each  group.  The  rows  of  these  matrices 
correspond  to  the  harmonic  number,  and  the  columns  correspond  to  each  (n,  m)  pair 
associated  with  that  group.  The  value  of  an  element  of  the  matrix  is  the  maximum 
value  of  the  C  parameters  of  all  the  data  points  Xj  contained  in  the  set  that  corresponds 
to  the  appropriate  group  and  (?z,  m)  pair.  This  matrix  shall  be  called  a  feature  matrix. 

5  Feature  Vectors  for  Classification 

The  feature  vector,  which  consists  of  a  set  of  pairs  of  TCF  Fourier  component  frequencies 
(cycle  frequencies  a)  and  magnitudes  (cyclic  cumulant  magnitudes  |CTCFQ|),  can  be  viewed 
as  a  one-dimensional  vector  with  complex-valued  elements  (a-f  i |CTCFa |)  indexed  by  (n,  m), 
or  as  a  two-dimensional  vector  (matrix)  with  real-valued  elements  indexed  by  [a,  (n,m)]. 
We  will,  therefore,  use  the  terms  feature  vector  and  feature  matrix  interchangeably.  The 
uniqueness  of  a  particular  signal’s  feature  vector  is  most  easily  appreciated  visually  by  using 
the  matrix  interpretation.  The  columns  of  the  feature  matrix  correspond  to  the  distinct 
order/number-of-conjugations  pairs  used  (the  ( n,m )  pairs),  and  the  rows  correspond  to  the 
harmonic  number  of  the  detected  harmonic  k,  where  a  —  k/To  for  which  T0_1  is  the  symbol 
rate  of  the  signal  that  gives  rise  to  the  feature.  The  value  of  the  matrix  for  a  given  row  and 
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column  is  the  maximum  (over  delay  sets)  absolute  value  of  the  amplitude  of  the  sine  wave 
with  the  harmonic  specified  by  the  row  and  corresponding  to  the  (n,m)  pair  specified  by 
the  column.  For  example,  the  GSA  might  produce  the  cycle  frequency  a0  for  order  2  and 
number  of  conjugations  1  for  several  different  delay  pairs  (t,-,  7y.).  The  GA  uses  all  of  these 
instances  of  detection  to  properly  group  all  the  GSA-produced  cycle  frequencies,  but  it  only 
outputs  the  largest  amplitude.  If  the  number  and  range  of  the  delay  vectors  that  are  input  to 
the  GSA  are  large  enough,  then  the  GSA  should  produce  among  these  instances  of  detection 
of  this  harmonic  one  that  corresponds  to  its  maximum  theoretical  value.  Ideally,  then,  the 
feature  matrix  consists  of  the  theoretical  maximum  of  the  amplitude  of  each  harmonic  that 
that  signal  actually  exhibits  for  each  input  (n,  m )  pair.  In  practice,  it  consists  of  an  estimate 
of  this  matrix.  The  definition  of  the  feature  matrix  is  shown  graphically  in  Figure  1. 


: 

: 

: 

*  ,  .  . 

k  =  2 

2/c  ±  2/To 

0  fc  ±  2/To 

k  =  1 

2 fc  ±  1/To 

0 fc  ±  1/To 

.  .  . 

k  =  0 

2/c 

0  fc 

(n,m)  =  (2,0)  (n,m)  =  (2,1) 


Figure  1:  Definition  of  a  feature  matrix  for  detection  and  classification.  An  element  of  the 
feature  matrix  corresponds  to  the  maximum  of  the  magnitude  of  the  cyclic  cumulant  for 
cycle  frequency  a  =  (n  —  2 m)fc  ±  k/T0,  where  k  is  uniquely  specified  by  the  row  and  the 
values  of  n  and  m  are  uniquely  specified  by  the  column. 

Examples  of  measured  feature  vectors  for  a  maximum  order  of  eight  (even  orders  only)  are 
presented  for  a  large  number  of  modulated  signals  of  interest  in  Figures  2-15.  These  examples 
include  PSK,  digital  QAM,  CPFSK,  and  partial-response  signaling.  These  features  were 
measured  from  data  records  of  length  16384  (1024  symbols)  samples  of  noise-free  simulated 
signals.  For  all  the  feature  matrices  shown  in  this  report,  the  first  three  columns  corresponds 
to  order  two,  the  next  five  correspond  to  order  four,  the  next  seven  to  order  six,  and  the 
remaining  nine  to  order  eight.  For  a  single  order,  the  columns  start  with  m  =  0  and  end 
with  m  =  n.  It  can  be  seen  from  these  examples  that  the  feature  vectors  are  distinct  for 
many  kinds  of  modulations,  but  that  a  sufficiently  large  order  of  processing  must  be  used 
to  insure  this  distinction.  For  example,  if  the  order  of  processing  is  restricted  to  six  or  less, 
then  the  feature  matrices  for  8PSK  and  16PSK  will  be  scaled  versions  of  each  other.  But 
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for  processing  order  of  eight  and  greater,  the  feature  matrices  for  these  two  modulations  are 
distinct. 

To  perform  classification  of  signals  using  these  features,  they  must  be  made  invariant 
to  the  signal  power  and  to  the  specific  carrier  frequency  and  symbol  rate.  By  design,  the 
features  are  already  invariant  to  the  latter  two  parameters:  the  same  feature  occurs  for  any 
carrier  and  symbol  rate  because  it  is  the  existence  and  relative  magnitudes  of  the  symbol-rate 
harmonics  that  define  the  feature.  However,  two  signals  that  are  identical  except  for  their 
power  levels  will  not  yield  the  same  feature.  Moreover,  the  two  resulting  features  are  not  just 
multiples  of  each  other.  The  effect  of  scaling  a  signal  on  the  classification  feature  is  explained 
in  the  next  section.  Following  that,  in  Section  6,  a  classification  algorithm  is  proposed,  along 
with  a  method  of  normalizing  the  measured  feature  for  the  class  of  rectangular-pulse  digital 
QAM  signals.  Further  research  is  required  to  generalize  the  method  to  other  signal  classes. 

5.1  Features  of  Scaled  Signals 

In  this  section,  the  effect  of  signal  scaling  on  the  classification  features  is  investigated.  This 
is  an  important  topic  because  the  power  of  the  signals  in  the  data  is  unknown.  The  effect  of 
scaling  a  signal  on  the  feature  matrix  is  not  linear.  In  particular,  the  scaling  of  the  signal  can 
be  interpreted  as  a  scaling  of  the  symbol  constellation  and,  therefore,  the  effect  of  scaling  is 
reflected  entirely  in  the  cumulant  of  the  symbol  variable.  Let  a  be  a  random  variable  with 
cumulants  C™n,  where  n  is  the  order  of  the  cumulant  and  m  specifies  the  number  of  times  the 
variable  is  conjugated  [11].  Then  the  random  variable  x  —  Ba ,  where  B  is  a  complex-valued 
constant,  has  cumulants  Bn~m{B*)mCfn. 

Let  s(t)  represent  a  complex-valued  PAM  signal  of  interest  with  symbol  sequence  {aj}. 
The  magnitudes  of  the  cyclic  cumulants  of  the  signal  x(t )  =  Bs(t ),  where  B  is  a  real  number, 
are  given  by 

Cm  r° o  n 

|G?(t)„|  =  \B\n  n  it  (2) 

to  1 

where  1  /T0  is  the  symbol  rate,  p{t )  is  the  pulse  function,  and  is  the  cycle  frequency  [11], 
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6  Classification  Algorithms 

6.1  Classification  From  the  GSA/GA  Output 

The  output  of  the  GSA/GA  combination  can  be  used  to  classify  the  corresponding  signal  in 
the  following  way: 

1.  Input  a  measured  feature  matrix  (one  of  the  groups  output  by  the  GA). 

2.  Compute  or  retrieve  theoretical  feature  matrices  for  all  desired  signal  types  for  the 
same  n,m  values  used  to  obtain  the  measured  feature  matrix. 

3.  Using  a  second-  or  fourth-order  element  of  the  measured  feature  matrix,  determine  the 
scaling  factor  (amplitude  of  the  signal).  (A  method  for  doing  this  for  rectangular-pulse 
digital  QAM  signals  is  explained  in  Section  6.2.)  Scale  the  measured  feature.  Ideally, 
the  candidate  matrix  for  the  true  signal  type  and  the  input  matrix  will  now  have  closely 
matching  element  values,  whereas  the  other  candidate  matrices  will  have  values  that 
are  significantly  different. 

4.  Compute  the  difference  between  the  normalized  input  feature  and  each  theoretical 
candidate  feature. 

5.  Declare  the  signal  to  be  of  the  type  that  produced  the  matrix  with  the  smallest  differ¬ 
ence  between  measured  and  theoretical  features. 

6.2  Rectangular-Pulse  Signaling 

In  this  section,  the  theoretical  feature  matrices  for  a  subclass  of  the  class  of  digital  QAM 
signals  are  computed  and  the  distance  between  them  is  evaluated.  The  subclass  is  that  for 
which  the  pulse  envelope  is  rectangular.  The  seven  signals  that  are  considered  are  BPSK, 
QPSK,  8PSK,  MPSK,  8QAM,  16QAM,  and  64QAM  (for  MPSK,  M  >  8).  These  signals 
contain  subsets  that  have  similar  cyclostationarity  properties  for  orders  1  through  8,  as  can 
be  seen  from  the  measurements  presented  in  Section  5.  The  PSK  signals  have  circular  symbol 
constellations  whereas  the  QAM  signals  do  not. 
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The  theoretical  feature  matrices  can  be  computed  for  this  subclass  because  the  maxi¬ 
mum  of  the  cyclic  temporal  cumulant  functions  is  relatively  easy  to  compute  for  all  orders 
and  cycle  frequencies.  In  particular,  for  any  full-duty-cycle  rectangular-pulse  PAM  signal 
with  independent  identically  distributed  symbols  {nm},  the  maxima  of  the  cyclic  temporal 
cumulant  functions  are  given  by 

max \CkJTa{T)n\  =  0  (3) 

which  occurs  whenever  the  difference  between  the  minimum  and  maximum  values  of  the 
elements  of  the  delay  vector  r  equals  To  —  (2 j  +  1)Tq/ (2|h|)  for  j  =  0,  •  •  • ,  \k\  —  1.  Note  that 
the  cumulant  C™n  incorporates  the  optional  conjugations  employed  in  Ck/T°(r)n.  For  k  =  0, 
the  maxima  are  equal  to  | C™n\  which  occur  whenever  the  delay  variables  are  equal  to  each 
other.  A  derivation  of  these  results  is  contained  in  Appendix  B.  Since  the  complex  envelope 
representations  of  digital  QAM  (including  PSK)  signals  are  exactly  this  type  of  PAM  signal 
(when  the  pulse  envelope  is  rectangular),  these  results  are  directly  applicable. 

The  result  (3)  can  be  used  to  construct  theoretical  feature  matrices  for  any  order  and 
number  of  harmonics.  In  this  section,  the  theoretical  features  for  the  above-mentioned  signals 
are  computed  and  the  differences  between  them  measured  for  the  purpose  of  evaluating  the 
potential  for  classification.  This  is  done  for  maximum  processing  orders  N  of  two  through 
eight.  The  results  are  displayed  in  terms  of  a  type  of  confusion  matrix.  This  confusion 
matrix  is  a  plot  of  the  distance  between  each  of  the  feature  matrices.  If  the  features  for  the 
various  signals  are  sufficiently  distinct,  then  the  confusion  matrix  will  be  zero  along  its  main 
diagonal  and  large  off  this  diagonal.  If  they  are  not  sufficiently  distinct,  then  small  values 
appear  off  the  main  diagonal.  The  confusion  matrices  are  shown  in  Figures  16-19.  These 
figures  show  that  the  feature  matrices  become  more  and  more  distinct  for  these  modulations 
as  the  order  increases. 

It  is  interesting  to  consider  the  possibility  of  performing  signal  classification  based  only 
on  the  symbol-rate  cyclic  features  (for  which  the  number  of  conjugated  factors  is  equal  to 
half  the  processing  order).  This  classification  possibility  is  well  motivated  because  the  most 
difficult  and  error-prone  part  of  measuring  the  feature  matrix  involves  the  estimation  of  the 
carrier-related  cyclic  features  (for  which  the  number  of  conjugated  factors  is  not  equal  to  half 
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the  processing  order).  The  difficulty  arises  from  the  fact  that  there  are  typically  many  more 
cycle  frequencies  for  n  ^  2m,  and  these  are  distinct  for  each  such  distinct  n  and  m.  There  are 
many  more  false  cycle  frequencies  that  appear  in  this  group  because  they  must  be  associated 
with  each  other  on  the  basis  of  their  distances  between  each  other.  By  coincidence,  there 
can  be  many  false  cycle  frequency  pairs  that  have  the  required  separation.  The  confusion 
matrices  for  the  reduced-dimension  classification  features  are  shown  in  Figures  20-23.  It  can 
be  seen  that  some  degree  of  classification  can  be  achieved  using  these  features,  but  that  it  is 
less  than  that  for  the  complete  features. 

6.3  Classification  From  Precise  Cumulant  Measurements 

The  GSA  and  GA  perform  the  task  of  estimating  and  grouping  the  cycle  frequencies  cor¬ 
responding  to  the  signals  in  the  data  quite  well,  but  the  GSA  is  prone  to  estimating  false 
cycle  frequencies.  These  false  cycle  frequencies  tend  to  degrade  the  quality  of  the  estimates 
of  higher-order  temporal  cumulant  functions,  which  manifests  itself  in  terms  of  errors  in  the 
estimate  of  the  strengths  of  the  cyclic  cumulants,  but  not  their  frequencies.  That  is,  the 
current  implementation  of  the  GSA  and  GA  algorithms  seems  to  be  most  well-suited  to  the 
cycle-frequency  estimation  problem,  not  the  cyclic-cumulant  estimation  problem. 

The  classifier  proposed  herein  depends  on  accurate  measurements  of  not  only  the  cycle 
frequencies,  but  also  the  cyclic  cumulants.  If  another  stage  of  processing  is  allowable,  the  cy¬ 
clic  cumulants  can  be  measured  using  previously  developed  software  (the  relevant  programs 
are  called  ctcfn,  partitions,  and  cross_cumulant,  which  are  written  in  C  and  were  de¬ 
veloped  by  the  first  author  during  his  doctoral  research  [4])  written  expressly  for  the  purpose 
of  estimating  cyclic  cumulants.  This  software  requires  knowledge  of  the  lower-order  cycle 
frequencies,  which  is  exactly  what  is  obtained  by  the  GSA/GA  tandem.  In  other  words,  the 
output  of  the  GSA/GA  tandem  can  be  used  by  ctcfn  and  cross_cumulant  to  perform  more 
precise  estimates  of  the  cyclic  cumulants,  which  will  result  in  more  accurate  estimates  of  the 
feature  vectors  and,  therefore,  better  classification. 
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7  Research  and  Development 

The  following  list  of  tasks  are  suggested  by  the  research  that  is  documented  in  this  report: 

1.  Computer  simulations  for  signal  environments  of  interest, 

2.  Mathematical  analysis  of  the  higher-order  cyclostationarity  of  communication  signals 
that  are  not  well  modeled  as  complex-valued  PAM  signals  at  baseband,  such  as  CPFSK, 
FSK,  and  analog  FM, 

3.  Cataloging  of  the  higher-order  cyclic  features  of  signals  of  interest, 

4.  Mathematical  characterization  of  the  quality  of  the  GSA’s  cyclic  cumulant  estimates 
as  a  function  of  SNR,  SIR,  and  collect  time, 

5.  Study  of  the  various  metrics  that  can  be  used  to  measure  the  distance  between  feature 
vectors  for  classification, 

6.  Development  of  generalizations  of  the  GA  to  handle  the  case  in  which  signals  share 
certain  cyclostationarity  properties  (e.g.,  two  or  more  signals  with  the  same  symbol 
rate  but  distinct  carrier  frequencies), 

7.  Design  and  development  of  fully  automated  and  operator-assisted  modes  of  operation, 

8.  Develop  alternative  hardware/software  architectures  for  the  two  modes  of  operation, 

9.  Specify  graphics  for  operator-assisted  mode  (and  for  demonstration), 

10.  Study  feasibility  of  integrating  into  the  GSA/GA  system  either  or  both  of  W.  A. 
Brown’s  MLSC  (maximum-likelihood  spectral  coherence)  classifier  and  W.  A.  Gard¬ 
ner’s  ASCR  (adaptive  spectral  correlation  recognition)  classifier,  both  of  which  utilize 
only  second-order  cyclostationarity, 

11.  Study  alternate  feature-extraction  scheme  that  uses  the  cycle  frequencies  computed  by 
the  GSA  and  grouped  by  the  GA  to  do  precise  cumulant  measurements  (see  Section 
6.3. 
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Figure  2:  Measured  feature  matrix  for  a 
rectangular-pulse  BPSK  signal. 
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Figure  3:  Measured  feature  matrix  for  a 
rectangular-pulse  QPSK  signal. 
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Figure  4:  Measured  feature  matrix  for  a 
rectangular-pulse  8PSK  signal.  This  feature 
contains  some  false  cycle  frequencies.  In  par¬ 
ticular,  the  only  upper  cycle  frequencies  for 
this  signal  and  n  <  8  consist  of  those  for  0 
and  8  conjugations. 
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Figure  5:  Measured  feature  matrix  for  a 
rectangular-pulse  offset  QPSK  signal. 


Figure  6:  Measured  feature  matrix  for  a  100%  Figure  8:  Measured  feature  matrix  for  a  25% 
excess-bandwidth  BPSK  signal.  excess-bandwidth  BPSI<  signal.  The  gray 

cells  for  n  =  4  and  6  and  large  k  are  false 
cycle  frequencies. 


Figure  7:  Measured  feature  matrix  for  a  50%  20 

excess-bandwidth  BPSK  signal.  The  gray 

cells  for  n  =  2  and  large  k  are  false  cycle  Figure  9:  Measured  feature  matrix  for  a 
frequencies.  duobinary  partial-response  PAM  signal. 


Figure  10:  Measured  feature  matrix  for  a 
CPFSK  signal  with  modulation  index  of  0.7. 
This  signal  does  not  fit  the  complex-valued 
PAM  model  that  the  GA  is  configured  to  rec¬ 
ognize,  and  so  the  feature  is  not  correct. 


Figure  12:  Measured  feature  matrix  for  a 
CPFSK  signal  with  modulation  index  of  1.0. 
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Figure  11:  Measured  feature  matrix  for  a 
CPFSK  signal  with  modulation  index  of  0.5. 
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Figure  13:  Measured  feature  matrix  for  a  4- 
QAM  signal.  Note  similarity  to  the  QPSK 
feature. 
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Figure  14:  Measured  feature  matrix  for  a  8- 
QAM  signal. 


Figure  15:  Measured  feature  matrix  for  a  16- 
QAM  signal. 
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Figure  16:  Theoretical-feature  confusion  ma¬ 
trix  for  seven  full-duty-cycle  rectangular- 
pulse  digital  QAM  signals  for  maximum  or¬ 
der  N  —  2  and  m  =  0,1,2.  The  signals 
are  (from  left  to  right  and  from  lower  to  up¬ 
per)  BPSK,  QPSK,  8PSK,  MPSK,  4QAM, 
8QAM,  16QAM,  and  64QAM. 


Figure  18:  Theoretical-feature  confusion  ma¬ 
trix  for  seven  full-duty-cycle  rectangular- 
pulse  digital  QAM  signals  for  maximum  or¬ 
der  N  —  6  and  m  =  0,  •  •  • ,  6.  The  signals 
are  (from  left  to  right  and  from  lower  to  up¬ 
per)  BPSK,  QPSK,  8PSK,  MPSK,  4QAM, 
8QAM,  16QAM,  and  64QAM. 


Figure  17:  Theoretical-feature  confusion  ma¬ 
trix  for  seven  full-duty-cycle  rectangular- 
pulse  digital  QAM  signals  for  maximum  or¬ 
der  N  —  4  and  m  =  0,  •  •  • ,  4.  The  signals 
are  (from  left  to  right  and  from  lower  to  up¬ 
per)  BPSK,  QPSK,  8PSK,  MPSK,  4QAM, 
8QAM,  16QAM,  and  64QAM. 


Figure  19:  Theoretical-feature  confusion  ma¬ 
trix  for  seven  full-duty-cycle  rectangular- 
pulse  digital  QAM  signals  for  maximum  or¬ 
der  N  =  8  and  m  =  0,  •  •  • ,  8.  The  signals 
are  (from  left  to  right  and  from  lower  to  up¬ 
per)  BPSK,  QPSK,  8PSK,  MPSK,  4QAM, 
8QAM,  16QAM,  and  64QAM. 
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Figure  20:  Theoretical-feature  confusion  ma¬ 
trix  for  seven  full-duty-cycle  rectangular- 
pulse  digital  QAM  signals  for  maximum  or¬ 
der  N  =  2  and  n  =  2m.  (The  features  consist 
only  of  symbol-rate  harmonics.)  The  signals 
are  (from  left  to  right  and  from  lower  to  up¬ 
per)  BPSK,  QPSK,  8PSK,  MPSK,  8QAM, 
16QAM,  and  64QAM. 


Figure  22:  Theoretical-feature  confusion  ma¬ 
trix  for  seven  full-duty-cycle  rectangular- 
pulse  digital  QAM  signals  for  maximum  or¬ 
der  N  =  6  and  n  =  2m.  (The  features  consist 
only  of  symbol-rate  harmonics.)  The  signals 
are  (from  left  to  right  and  from  lower  to  up¬ 
per)  BPSK,  QPSK,  8PSK,  MPSK,  8QAM, 
16QAM,  and  64QAM. 


Figure  21:  Theoretical-feature  confusion  ma¬ 
trix  for  seven  full-duty-cycle  rectangular- 
pulse  digital  QAM  signals  for  maximum  or¬ 
der  N  =  4  and  n  —  2m.  (The  features  consist 
only  of  symbol-rate  harmonics.)  The  signals 
are  (from  left  to  right  and  from  lower  to  up¬ 
per)  BPSK,  QPSK,  8PSK,  MPSK,  8QAM, 
16QAM,  and  64QAM. 


Figure  23:  Theoretical-feature  confusion  ma¬ 
trix  for  seven  full-duty-cycle  rectangular- 
pulse  digital  QAM  signals  for  maximum  or¬ 
der  N  =  8  and  n  =  2m.  (The  features  consist 
only  of  symbol-rate  harmonics.)  The  signals 
are  (from  left  to  right  and  from  lower  to  up¬ 
per)  BPSK,  QPSK,  8PSK,  MPSK,  8QAM, 
16QAM,  and  64QAM. 
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A  The  Theory  of  Higher-Order  Cyclostationarity 

In  this  appendix  we  define  the  temporal  and  spectral  moment  and  cumulant  functions  that 
form  the  basis  of  the  theory  of  HOCS.  Then  we  give  a  brief  tutorial  explanation  of  how  cum- 
ulants  arise  as  the  solution  to  the  problem  of  generating  pure  nth-order  sine  waves.  Finally, 
we  explain  the  signal-selectivity  property  that  is  unique  to  the  cyclic  temporal  cumulants 
and  their  Fourier  transforms,  the  cyclic  polyspectra,  and  we  illustrate  the  parameters  using 
the  example  of  digital  QAM  signals. 

For  a  time-series  x{t)  for  — oo  <  t  <  oo,  we  define  the  nth-order  lag-product  time-series 

by 

Lx{t,r)n=  f[  x(t  +  Tj),  (4) 

j=i 

where  t  =  [rj  •  •  •  rn]f  and  [•]!  denotes  matrix  transposition.  The  cyclic  temporal  moment 
function  (CTMF)  of  order  n  is  defined  by  the  limiting  time  average 

K(r)n  =  Hm  =  (5) 
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which  is  simply  the  Fourier  coefficient  associated  with  the  component  el2lzat  in  the  time- 
series  Lx(t,r)n.  It  can  be  seen  that  the  CTMF  is  a  Fourier  coefficient  of  a  moment  function 
because  the  nth-order  fraction-of-time  probabilistic  moment  (the  temporal  moment  function 
[TMF])  associated  with  the  lag  product  Lx(t,r)n  is  given  by 

R,(t,r)n=  £(“)  T)„}  (6) 

where  {•}  =  ({-)e~l2*at}  e,2™‘,  and  can  be  expressed  as  [14,  15] 

R,(t,T)„=  E  (7) 

og{cr} 

where  the  sum  is  over  all  real  numbers  a ,  called  nth-order  cycle  frequencies ,  for  which 
Rx(r)n  ^  0.  In  (6),  E{a}{-}  is  the  temporal  expectation  operation  (or  the  sine-wave  extrac¬ 
tion  operation).  The  functions  (5)  and  (6)  exist  and  are  well-behaved  for  appropriate  mod¬ 
els  of  many  time-series  including  amplitude  modulated  (AM),  pulse-amplitude-modulated 
(PAM),  phase-shift-keyed  modulated  (PSK),  and  digital  quadrature  AM  (QAM)  signals 
[10,  4],  and  others. 

The  temporal  cumulant  function  (TCF)  of  order  n  for  the  set  of  time-series  translates 
{x(i  +  Tj)}”=1  is  defined  by 

Cx(t,T)n  =  ^  k(p)  Rx{ti  TvJ  )n:  5  (8) 

P=M Li  L 

which  is  completely  analogous  to  its  stochastic-process  counterpart  in  [16].  The  sum  in  (8)  is 
over  all  distinct  partitions  Pn  of  the  set  of  indices  {1, 2,  •  •  ■ ,  n},  where  each  partition  {uk}vk=i 
has  p  elements,  1  <  p  <  n,  and  k(p)  =  (-l)p-1(p  -  1)!.  The  vector  rU]  is  the  vector  of  nj 
lags  with  indices  in  the  set  Uj. 

The  cyclic  temporal  cumulant  function  (CTCF)  of  order  n  is  the  Fourier  coefficient  of  the 
TCF: 

Cf(r)„=  r)ne-2^‘).  (9) 

The  set  of  real  numbers  {/?}  for  which  Cf(r)„  ^  0  is  called  the  set  of  pure  nth-order  cycle 
frequencies ,  for  reasons  that  will  become  clear  subsequently.  Combining  (6)-(9)  reveals  that 
the  CTCF  is  given  by  the  following  explicit  function  of  lower-order  CTMFs: 

cf(r)„  =  E  kp)  E  (n«?KM  -  <10) 

p  L  at  1=(3  \1=1  /. 

where  1  =  [1  •  •  •  1]+,  and  a  =  [ai-'-cq,]'1'.  The  CTCF  was  originally  derived  in  [2]  as  the 
solution  to  the  problem  of  removing  from  the  Fourier  coefficient  all  contributions 

from  Fourier  coefficients  Rx3{rUj)nj  of  lower  order.  This  is  equivalent  to  removing  from  the 
finite-strength  additive  sine-wave  component  of  frequency  (3  in  the  lag  product  time-series 
Lx(t ,  r)n  all  contributions  from  products  of  sine-wave  components  in  lag  products  Lx{t ,  tVj  )nj 
of  lower  order — whose  frequencies  sum  to  (3 — that  can  be  obtained  by  factoring  Lx(t,T)n. 
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By  using  the  relationship  between  moments  and  cumulants[21],  the  CTMF  can  be  expressed 
in  terms  of  CTCFs  of  order  n  and  lower: 


R°r(T)n 


E 

P 


e  n^K), 

lp'l=aj=l 


(11) 


The  CTMF  and  the  CTCF  are  not  in  general  integrable  due  to  the  presence  of  sinusoidal 
components  in  r.  These  components  formally  result  in  Dirac  deltas  in  the  n-dimensional  Fou¬ 
rier  transform  of  the  CTCF.  However,  a  reduced-dimension  version  of  the  CTCF  is  absolutely 
integrable  for  many  time-series  of  interest  and,  therefore,  it  is  strictly  Fourier  transformable 
[2].  The  reduced-dimension  CTCF  (RD-CTCF)  is  simply  the  CTCF  associated  with  the  n 
variables  {x(t  +  Tj)}”=i  with  rn  =  0.  The  RD-CTCF  is  denoted  by 

££(“)«=  Cx(T)n  for  n  =  Ui  and  rn  =  0,  (12) 


where  wis  (n-l)-dimensional  vector  [ui  i].  The  (n  -l)-dimensional  Fourier  transform 

of  (12)  is  denoted  by  P£(f')n: 


P,0(f)n  =  r  ■■■  P  Cl («)„ e-a”u'f  in,  (13) 

J — oo  J — oo 


where  /'  =  [fi  •  •  •  i]f. 

The  spectral  moment,  spectral  cumulant,  and  cyclic  polyspectrum  are  defined  as  follows. 
Consider  the  n  complex-demodulate  time-series  )  for  j  =  1  associated  with 

narrow  bandpass  filtered  versions  of  x(t),  where 

rt+T/2 

XT(tJ)=  x(v)e-'2”fvdv.  (14) 

Jt-T/2 

The  limit  as  T  — >  oo  of  the  limiting  time-average  of  the  product  of  these  spectral  components 
is  called  the  spectral  moment  function  (SMF)  of  order  n 

&(/)»=  Tlto  (f[Vx(i, /■>)),  (15) 

and  it  can  be  shown  that  Dirac  deltas  can  be  factored  out  as  follows: 

Si(/)„  =  E5,°(/%'S(/,1-q).  (16) 

a 

where  S(-)  is  the  Dirac  delta  function.  However,  the  factor  Sf(f)n  contains  additional  Dirac 
deltas  for  many  signals  and  n  >  2  (e.g.,  BPSK  and  n  —  4). 

The  spectral  cumulant  function  (SCF)  of  order  n  is  given  by 


Px{f)n 


E 


twns.(/,,). 

j=i 


(17) 
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where  fu.  is  the  vector  of  nj  frequencies  with  subscripts  in  the  set  i /j,  and  it  follows  from 
(16)  that  Dirac  deltas  can  again  be  factored  out: 

ft(/)»  =  EWM(/t  1-P)-  (is) 

ft 

Analogous  to  the  definition  of  the  cyclic  spectrum  (and  power  spectrum)  in  [14],  the 
factor  Pf(f')n  is  defined  to  be  the  cyclic  polyspectrum  (CP)  and  is  given  explicitly  by 


£?(/')»  =  £ 


(19) 


As  first  shown  in  [2],  the  CP  is  the  (n  -  l)-dimensional  Fourier  transform  of  the  RD-CTCF 
C^(u)n  (cf.  (13)).  This  is  a  generalization  of  the  Wiener  relation  between  the  power  spec¬ 
trum  and  autocorrelation  from  second-order  stationary  time-series  (cf.  [14])  to  nth-order 
cyclostationary  time-series.  Within  the  stochastic-process  framework  of  generally  nonsta¬ 
tionary  processes,  it  should  be  called  the  cyclic  Shiryaev-Kolmogorov  relation  [16],  which  is 
the  generalization  of  the  Wiener-Khinchin  relation  (cf.  [14]).  Because  the  RD-CTCF  can  be 
absolutely  integrable,  the  CP — unlike  the  SMF — contains  no  Dirac  deltas.  That  is,  all  Dirac 
deltas  present  in  the  individual  terms  in  (19)  cancel. 

Let  us  now  see  how  the  cumulant  arises  as  the  solution  to  the  problem  of  generating 
pure  nth-order  sine  waves.  For  low  orders  n,  it  is  easy  to  mathematically  characterize  a 
pure  nth-order  sine  wave  in  a  way  that  matches  our  intuition.  For  notational  simplicity, 
we  choose  the  case  of  no  conjugated  factors  in  (4).  For  n  =  1,  the  moment  sine  waves 
(the  sine-wave  components  of  the  polyperiodic  TMF)  are,  by  definition,  pure  lst-order  sine 
waves.  For  n  =  2,  all  products  of  lst-order  moment  sine  waves  can  be  subtracted  from  the 
2nd-order  moment  sine  waves  to  obtain  the  pure  2nd-order  sine  waves,  which  are  denoted 
by  ax(t,Ti,  r2)2  , 


crx(t,Ti,T2)2  =  E{a}  {x(t  +  Ti)x(t  +  r2)}  -  E{a}  {x(t  +  n)}  £{a}  {x(t  +  r2)} 
=  Rx(Et)2  -  Rx{t,Tl)iRx{t,T2)1. 


There  are  several  interesting  points  to  be  made  concerning  pure  2nd-order  sine  waves: 
( i )  since  Rx(t,  r,)i,  i  =  1,2,  and  Rx(t,r) 2  are  1st-  and  2nd-order  moments  respectively,  then 
ax(t,  Ti,  t2) 2  is  a  temporal  covariance  function;  ( ii )  if  Rx(t ,  r)i  =  0,  then  there  are  no  lower- 
order  sine  waves,  and  the  2nd-order  moment  sine  waves  are  equal  to  the  pure  2nd-order 
sine  waves;  (Hi)  if  the  variables  x(t  +  tj)  and  x(t  +  t2)  are  statistically  independent  (in  the 
temporal  sense  [14,  15]),  then  E^  {x(t  +  T])x(t  +  r2)}  =  E^  {x(t  +  Ti)}  E ^  {x(t  +  r2)} 
and  therefore  ax(t ,  rl5  r2)2  =  0,  that  is,  there  is  no  pure  2nd-order  sine  wave  for  this  particular 
pair  of  delays  and  r2.  This  latter  point  shows  that  the  set  of  pure  cycle  frequencies  {/3} 
can  be  considerably  different  from  the  set  of  impure  (moment)  cycle  frequencies  {a}. 

The  pure  3rd-order  sine  waves  are  obtained  next.  From  the  3rd-order  moment  sine  waves, 
we  want  to  subtract  each  possible  product  of  lower-order  sine  waves,  but  only  once  each. 
Thus,  we  subtract  products  of  pure  2nd-order  and  pure  lst-order  sine  waves  from  the  3rd- 
order  moment  sine  waves,  rather  than  subtracting  products  of  1st-  and  2nd-order  moment 
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sine  waves: 


crx(t,r) 3  =  £{a}  x(f  +  Tj)|  -  <Ja;(^r1,T2)2  crx(i,r3)i  -  crx{t,  tu  r3)2  ax(t,  t2)i 

0jf(^j  T2,  ^"3)2  @x(tl  P  )  1  ^x(^,^"l)l  ^xij-1  p)l  ^3)1  ■ 

Observe  that  there  are  no  other  possible  products  of  pure  lower-order  sine  waves.  The  terms 
in  the  sum  of  products  that  are  subtracted  can  be  enumerated  by  considering  the  distinct 
partitions  of  the  index  set  {1,2,3}.  A  partition  of  a  set  G  is  a  collection  of  p  subsets  of  G, 
{vjYj- j,  with  the  following  properties:  G  =  (Jj=i  and  vj  f|  *4-  =  0  for  j  /  k.  The  set  P3 
of  distinct  partitions  of  {1,2,3}  is  given  by 

P  =  1  :  {1,2,3} 

P  =  2:  {1,2},  {3}  {1,3},  {2}  {2,3},  {1} 

P  =  3:  {1},  {2},  {3}. 

Thus,  we  can  express  the  pure  3rd-order  sine  wave  <7x(f,r)3  as  a  sum  over  the  elements  of 
Pz: 


<7x(t,r)3  =  Rx(t,T)3 


p3 

P^l 


13  <rx(t,  tUj  )Uj 

3= 1 


(20) 


where  tUj  is  a  subset  of  {tj}|=1  with  elements  having  indices  in  vj,  and  rij  is  the  number  of 
elements  in  Uj.  Notice  that,  as  in  the  case  of  n  =  2,  if  the  lst-order  moment  sine  waves  are 
zero,  then  the  3rd-order  moment  sine  waves  are  equal  to  the  pure  3rd-order  sine  waves.  In 
this  case,  there  are  no  products  of  lower-order  sine  waves  to  subtract  from  the  moment  sine 
waves. 

The  formula  for  the  pure  nth-order  sine  waves  is 


(Tx(t,T)n  =  Rx{t,r)n  -J2 


Pn 
P#  1 


31  {t->  Tv}  )rij 

i=i 


(21) 


where  Pn  is  the  set  of  distinct  partitions  of  the  index  set  {1,2,  •  •  ■  ,n}.  The  pure-sine-waves 
formula  (21)  gives  all  the  pure  nth-order  sine  waves  associated  with  the  delay  set  r.  A 
single  pure  ?zth-order  sine  wave  with  frequency  (3  can  be  selected  by  computing  the  Fourier 
coefficient: 

^(r)„  ei2^  =  (ax(u,  r)B  e~i2^u^)  ,  (22) 

and  can  be  expressed  in  terms  of  pure  lower-order  sine  waves  by  using  the  Fourier  series  for 
each  <Jx{t,TVj)nj  in  (21): 


ax(t,w)k  =YJ°xk(w)ket'27r0kt,  w  =  [wi  ■  ■  ■  tu*]!, 

Pk 


(23) 


where  the  sum  is  over  all  cycle  frequencies  (3k  of  order  k.  Thus,  the  strength  of  a  single  pure 
nth-order  sine  wave  is  given  by 


°x(T)n  =  RPx{r)n  -  J2 


Pn 
P#  1 


J2  II  ax(T^)n3 


ti  LI 'P=p  i_1 


(24) 
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where  (3  is  the  p-dimensional  vector  of  pure  cycle  frequencies  [(3\  •  ■  ■  /3p]*  and  1  is  the 
p-dimensional  vector  of  ones.  Hence,  the  pure-sine-wave  strength  cr^(r)n  is  given  by  the 
CTMF  RP(r)n  with  all  products  of  pure  lower-order  sine-wave  strengths,  for  sine  waves 
whose  frequencies  sum  to  /?,  subtracted  out.  From  (11),  it  is  clear  that  the  pure-sine-wave 
strength  cr^(r)n  is  identical  to  the  CTCF  Cf(r)n. 

Finally,  let  us  see  how  the  CTCF  (and  the  CP)  is  signal  selective.  Let  the  signal  x(t) 
consist  of  the  sum  of  M  mutually  independent  signals  u 

M 

z(t)  =  y ■»»(<)■  (25) 

m= 1 

Then  the  addition  rule  for  cumulants  [4,  17]  can  be  used  to  show  that  the  nth-order  TCF 
for  x(t)  is  the  sum  of  TCFs  for  {ym(2)}> 

M 

Cx(t,r)n  =  Cym(t,r)n.  (26) 

771  =  1 

Thus,  the  pure  nth-order  sine  waves  in  the  delay  products  of  each  ym(t)  add  to  form  the 
pure  nth-order  sine  wave  in  the  delay  product  of  x(t): 

M 

Cf(T)„  =  •£  Cl(T)„.  (27) 

m=l 

The  TMF  does  not  in  general  obey  this  useful  cumulative  relation  (27).  That  is,  the  nth-order 
TMF  for  x(t)  is  not  the  sum  of  nth-order  TMFs  for  each  ym(t ):  Rx(t,  r)n  ^  Y2m=i  Rym(ti  T)n 
.  An  exception  is  the  case  of  zero-mean  signals  and  n  <  3,  for  which  moments  and  cumulants 
are  equal,  which  is  a  commonly  encountered  case  in  HOS. 

To  illustrate  how  (27)  can  be  applied  in  practice,  consider  the  situation  in  which  {ym(t)} 
represent  M  interfering  signals  that  overlap  in  time  and  frequency,  but  each  ym(t)  possesses 
some  distinct  nth-order  cycle  frequency,  say  (3m.  Then  it  follows  from  (27)  that 

Cf”(r)„  =  C£(  r)„,  m  =  1, 2,  •  •  • ,  M.  (28) 

This  indicates  that  the  presence  or  absence  of  each  of  the  signals  ym(t)  can  be  detected 
by  measuring  (estimating)  the  CTCFs  of  x(t)  for  the  cycle  frequencies  {Rm},  and  that 
parameters  of  each  of  the  signals,  on  which  these  CTCFs  depend,  can  be  estimated.  As 
illustrated  in  [14,  12,  13]  for  second  order  and  in  [4]  for  higher  order,  this  signal-selectivity 
property  can  be  exploited  in  numerous  ways  to  accomplish  noise-and-interference-tolerant 
signal  detection  and  estimation. 

As  another  application,  let  M  =  2,  yi(t)  be  non-Gaussian,  and  j/2^ )  be  Gaussian.  Then 
Cy2(t,r)n  =  0  for  n  >  3  and,  from  (26),  we  have  Cx{t,r)n  =  C'yi (T r)„,  n  >  3,  which 
indicates  the  detectability  of  y\(t)  with  no  knowledge  about  y2(i)  except  that  it  is  Gaussian. 

As  an  example  of  the  cyclic  polyspectrum,  we  consider  the  class  of  digital  QAM  signals 
described  by  their  complex  envelopes,  which  can  be  expressed  as  x(t)  =  Z)m=-oo  amP{t  — 
niTo  —  to),  where  p(t)  is  the  complex  pulse  and  {am}  is  the  complex  digital  data.  Assuming 
that  { am }  is  an  independent  sequence,  we  have  shown  [10,  5,  11]  that  the  CP  for  the  set  of 
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variables  {z(*b(i  +  Tj)}"=1  (with  T„ 
jth  variable,  is  given  by 

P?(f%  =  %^r((-)n[/J  - 
-^0 


=  0),  where  (*)j  denotes  an  optional  conjugation  of  the 

■  lt/'])(*)B  n  ^((-)J/J)(*)jei2^°,  (3  =  k/T0,  (29) 

i=i 


where  (-)j  is  the  optional  minus  sign  corresponding  to  the  optional  conjugation  (*)j,  Ca)„ 
is  the  cumulant  of  am,  and 

/OO 

p(t)e-2’"  it.  (30) 

-oo 

Thus,  the  CP  is  simply  a  scaled  product  of  n  pulse  transforms. 


B  Cumulants  of  Rectangular-Pulse  PAM  Signals 


The  maxima  of  the  absolute  value  of  the  nth-order  cyclic  cumulants  for  PAM  signals  with 
full-duty-cycle  rectangular  pulses  are  derived  in  this  appendix. 

For  the  signal  model 

OO 

x(t)=  Y  aiP(t  +  jT0  +  t0),  (31) 


j=-oo 


where 


p(t)  = 


1,  |f|<7b/2 
0,  1<|  >  To/2, 


(32) 


and  {aj}  is  an  independent,  identically  distributed  symbol  sequence,  the  nth-order  cyclic 
temporal  cumulant  function  (CTCF)  is  given  by 

cm  r°°  n 

Cx(T )n  =  ^  n  p(<  +  Tk)e~i2*at  dt  ei2vat° .  (33) 

Jo  J-oo  k=1 

Let  t\  equal  the  difference  between  the  largest  and  smallest  of  the  delay  variables 

ti  =  max  | Ti  —  Tj  | 

if  this  number  is  less  than  T0  and  equal  to  T0  otherwise.  Then  the  product  of  the  shifted 
rectangle  functions  in  the  integrand  of  (33)  is  itself  equal  to  a  rectangle  with  width  T0  —  t\. 
It  is  easy  to  evaluate  the  integral  (33),  differentiate  the  result  with  respect  to  ti,  set  this 
derivative  equal  to  zero  and  solve  for  t\.  The  result  is 


max|C“(r)7 


|CQmn|  r  4  (2|fc|-2m-l)T0  „  „ 

for  U  =  - — - - — ,  m  =  0,  \k  —  1|. 


ir\k\ 


2\k\ 
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